{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概率判别式模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "三种分类方法：\n",
    "\n",
    "(1) Discriminant function（判别函数）：直接把输入 x 判归某一个类别\n",
    "\n",
    "(2) Generative model（生成模型）：通过计算 joint distribution 或 class-conditional distribution$p(x|C_k)$而计算 posterior $p(C_k|x)$\n",
    "\n",
    "(3)Discriminant model（判别模型）：直接计算 posterior $p(C_k|x)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3.2 Logistic 回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在⼀些相当⼀般的假设条件下，类别C1的后验概率可写成作⽤在特征向量$\\phi$的线性函数上的logistic sigmoid函数的形式:\n",
    "$$p(C_1|\\phi)=\\sigma(w^T\\phi)$$\n",
    "对于⼀个M维特征空间ϕ，这个模型有M个可调节参数。不同于最大似然方法确定高斯类条件概率密度时参数随M二次增长($\\frac{M(M+5)}{2}$)，这里对参数数量M呈线性依赖\n",
    "$$\\frac{d\\sigma}{da}=\\sigma(1-\\sigma)$$\n",
    "$$p(t|w)=\\prod_{n=1}^Ny_n^{t_n}(1-y_n)^{1-t_n}$$\n",
    "对似然函数去负对数就产生了交叉熵误差函数：\n",
    "$$E(w)=-ln(p|t)=-\\sum_{n=1}^N{t_nlny_n+(1-t_n)ln(1-y_n)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3.3 牛顿迭代法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最小化$E(w)$时令其梯度为0：$\\nabla E(w)=\\sum_{n=1}^N(y_n-t_n)\\phi_n=0$\n",
    "\n",
    "在第3章讨论线性回归模型的时候，在⾼斯噪声模型的假设的情况下，最⼤似然解有解析解。这是因为对数似然函数为参数向量w的⼆次函数。对于logistic回归来说，不再有解析解了，因为logistic sigmoid函数是⼀个⾮线性函数。然⽽，函数形式不是⼆次函数并不是本质的原因。精确地说，正如我们将要看到的那样，误差函数是凸函数，因此有⼀个唯⼀的最⼩值。此外，误差函数可以通过⼀种⾼效的迭代⽅法求出最⼩值，这种迭代⽅法基于Newton-Raphson迭代最优化框架，使⽤了对数似然函数的局部⼆次近似。\n",
    "\n",
    "$$w_{new}=w_{old}-H_{-1}\\nabla E(w)$$\n",
    "$$H=\\nabla\\nabla E(w)=\\sum_{n=1}^Ny_n(1-y_n)\\phi_n\\phi_n^T=\\Phi^TR\\Phi$$\n",
    "\n",
    "R为nxn的对角矩阵，其元素为$R_nn=y_n(1-y_n)$,这样更新过程变为：\n",
    "\n",
    "$$\\begin {align}\n",
    "{new}&=w_{old}-(\\Phi^TR\\Phi)^{-1}\\Phi^T(y-t)\\\\\n",
    "     &=(\\Phi^TR\\Phi)^{-1}\\Phi^TR\\Phi w_{old}-\\Phi^T(y-t) \\end {align}$$\n",
    "\n",
    "由于权矩阵R不是常\n",
    "量，⽽是依赖于参数向量w，因此我们必须迭代地应⽤规范⽅程，每次使⽤新的权向量w计\n",
    "算⼀个修正的权矩阵R。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(r\"../\")\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "from prml.features import PolynomialFeatures\n",
    "from prml.linear import (\n",
    "    BayesianLogisticRegressor,\n",
    "    LeastSquaresClassifier,\n",
    "    LinearDiscriminantAnalyzer,\n",
    "    LogisticRegressor,\n",
    "    Perceptron,\n",
    "    SoftmaxRegressor\n",
    ")\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def create_toy_data(add_outliers=False, add_class=False):\n",
    "    x0 = np.random.normal(size=50).reshape(-1, 2) - 1\n",
    "    x1 = np.random.normal(size=50).reshape(-1, 2) + 1.\n",
    "    if add_outliers:\n",
    "        x_1 = np.random.normal(size=10).reshape(-1, 2) + np.array([5., 10.])\n",
    "        return np.concatenate([x0, x1, x_1]), np.concatenate([np.zeros(25), np.ones(30)]).astype(np.int)\n",
    "    if add_class:\n",
    "        x2 = np.random.normal(size=50).reshape(-1, 2) + 3.\n",
    "        return np.concatenate([x0, x1, x2]), np.concatenate([np.zeros(25), np.ones(25), 2 + np.zeros(25)]).astype(np.int)\n",
    "    return np.concatenate([x0, x1]), np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYnGXV+PHvmb4lyWazKaQAyRLEQAglhCYQQjEJTZQu\nqKDygtgLUuz68r6W3yuoKEaUIiCi0gkl9A4ptAQCaaST3rZMP78/nkkys/vM7uzO7E7Z87muvdx5\n6g0uZ+5y7vsWVcUYYyqJp9gFMMaYQrPAZoypOBbYjDEVxwKbMabiWGAzxlQcC2zGmIpjgc0Y0+NE\n5G8isl5E5mc5LyLyOxFZLCJvi8ghaeemisj7qXNX5fI+C2zGmN5wKzC1g/PTgLGpn0uBPwGIiBe4\nMXV+HHC+iIzr7GUW2IwxPU5Vnwc2d3DJGcDt6ngVqBORPYBJwGJVXaqqUeDu1LUd8hWi0F1VP7Be\nR4wYVYxXF0UkHscf2P2vOpKME/Rn/qtXDRP0+nu7aKZCvTlvwUZVHZzPM044+RO6edPWXN+3AAin\nHZqhqjO68LoRwMq0z6tSx9yOH97Zw4oS2EaMGMW998wsxquLYsnGzYwaPSTzWPMGGofW7/ocji8C\nYEz/vP4WjQGgPrT/8nyfsXnTVp5++Z5c3xdW1Yn5vrNQrCnaS1YuW9/h+ZBvbC+VxJiStBpIb8aN\nTB3LdrxDFth6QWNDvevxJes66nIwpk95EPhcanT0CGCbqq4FZgNjRWS0iASA81LXdqgoTVEDjTWD\nWdK8IeNYyDeWpdsXWXPUVBwR+QcwGWgQkVXAjwE/gKreBMwEpgOLgRbg4tS5uIh8FXgc8AJ/U9UF\nnb3PAlsvaWyoZ8my9e372tZtzuhrM6YSqer5nZxX4Ios52biBL6cWVO0iBpr3GtmS7dvcD1ujMmN\nBbZeZoMIxvQ8C2y9yG0QobFmsOsggtXajOk+C2xFYLU2Y3qWBbZeZqkfxvQ8C2wlwG0QwUn9sOao\nMd1hga0IGhvqXZujVmszpjAssJUIS/0wpnAssBWRDSIY0zMssBWJpX4Y03MssBWZ1dqMKTwLbEVk\nqR/G9AwLbCXGUj+MyV/BApuIeEXkDRF5uFDP7Ass9cOYwitkje0bwHsFfF6fZakfxuSnIIFNREYC\npwA3F+J5fZENIhhTOIWqsV0PXAkks10gIpeKyBwRmbN586YCvbYyWOqHMYWVd2ATkVOB9ao6t6Pr\nVHWGqk5U1Yn19YPyfW1FslqbMYVRiBrb0cDpIvIhzmamU0TkjgI8t0+x1A9jCifvwKaqV6vqSFXd\nG2cHmadV9cK8S2Ys9cOYbrI8thJiqR/GFEZBA5uqPquqpxbymX2dpX4Y03VWYytBNohgTH4ssJUY\nS/0wJn8W2EqU1dqM6T4LbCXIUj+MyY8FtjJhqR/G5M4CW4my1A9jus8CWxmx1A9jcmOBrcR1Vmuz\nQQRj2rPAVsKypX64sVqbMbtZYCsDlvphTNdYYCtxXUn9sFqbMQ4LbGUivdaWLfXDGOOwwFYGip2w\nq5HX0K1Xo1suR5tuQhO2ArIpbRbYylRvDSJoy33QdD3E34fEOgg/BVu/jSa2FPQ9xhSSBbYy0tup\nH5pshZZ/gkbSjiZAW6H1/oK+y5hCssBWJoqS+pFYBeJ1ORGH+NuFeYcxPcACW5np1dQPTx1o3P2c\nuAdVY0qBBbYy0tupH+IdDP5xgK/NmSBUfyrv5xvTUyywlaFeTf3o913wjwd8ICGQaqi9FPGPK9w7\njCmwtl/FpsQ1NtSzZGP7GtqSdZtpHOpeo8uHeGpgwA/RxDbQ7eDdAxH7szGlzWpsFaA3Uj/EOwDx\njbKgZsqCBbYylMtabTYTwfRlFtgqRCWv+qGJbWjz7ejWb6Hbfo5G3yl2kUyJs8BWxvrCqh+a2Arb\nvgWtD0N8OcTegO3Xoa2PF7topoRZYCtTfWabvtb7IdkEpOfTRaDlNlRjxSqV6QYRmSoi74vIYhG5\nyuX890TkzdTPfBFJiEh96tyHIvJO6tyczt5lga3MVVKtTTWGRl5Bw4+i8WXOwdgbZAa1nQQSK3qz\neCYPIuIFbgSmAeOA80UkI2dIVX+tqgep6kHA1cBzqpr+TX186vzEzt5nQ1xlrLdTP3qSxlfCth8C\nMdCEcyxwKEgdsNLtBpABvVrGviaSiBWytj8JWKyqSwFE5G7gDODdLNefD/yjuy+zwFZhGmsGs6Q5\n84/R2aZvEWP6l+Y0KFWFHb9y8uTSRedB6ASIfwCkT8T3gm9fxNvQm8Xsc0RCXanxN7RpIs5Q1Rlp\nn0eQ+Q21Cjjc/b1SDUwFvpp2WIEnRSQB/LnNs9uxwFbmGhvqWbJsPaNGD8k4Xla1tsRHkHCrGUQg\nvhBqLoSWOwCPU1Pz7QP9r+ztUpqObcyliZij04CX2jRDP6Gqq0VkCDBLRBaq6vPZHmCBrQK51drA\nGUQozVpbHESc7+S2NI5UnYKGToT4CvDUOXNYTblZDYxK+zwydczNebRphqrq6tT/rheR+3CatlkD\nmw0eVIiyHkTwjnTmoLYTgOCxAIgEEf9YC2rlazYwVkRGi0gAJ3g92PYiERkAHAc8kHasRkT67fwd\nOBmY39HLLLBVgHJP/RAR6PcdZ5I9/tTBEPhGQdUpRS2bKQxVjeP0mT0OvAfco6oLROQyEbks7dIz\ngSdUtTnt2FDgRRF5C3gdeERVH+vofdYUrSArXfra0oV8YwnHF/ViiXIn/nFo3Y0QeRaSG8F/AAQm\nIa4LXZpypKozgZltjt3U5vOtwK1tji0FJnTlXRbYKkQlpH6IdyBUn1nsYpgKkHdTVERGicgzIvKu\niCwQkW8UomAmf9nWaivF5mhXaXwl2vx3tOmvaHS+kzJiTEohamxx4DuqOi/VwTdXRGaparbEO9ND\nKiL1Iwfa+ig034bzp5eEyJMQOBKt/ZrTX2f6vLxrbKq6VlXnpX7fgdMxOCLf55rCqLRVPzSxFZpv\nBaJAMnUwAtFXIGarfhhHQUdFRWRv4GDgtUI+13RNWad+dCb2hvvOWRqB6Mu9Xx5TkgoW2ESkFvgP\n8E3VtnNjQEQuFZE5IjJn82bbSbynlHvqR+f8gFtzU4BAL5fFlKqCBDYR8eMEtTtV9V63a1R1hqpO\nVNWJ9fWDCvFa04GKrbUFDmVXEzSDH0LH93ZpTIkqxKioAH8F3lPV/8u/SCZfXdmmr1RoYgPa+hDa\n+iAa/yjrdeKpgn5XAkGQqt1JvdXnIr7RvVZeU9oKMSp6NHAR8I6IvJk6dk0qGc+UiFJe9cMZ5byV\n3ZNF70RrLkCqznC9XgIHo/V/g+hsIAr+gxGvtQLMbnkHNlV9EfdOD1NE5ZL6oYkNqaDWZjXc5rtQ\n/2GIb7jrfeKpgtCxPV4+U55srmgfUpKpH9HXs5xIOikcxnSDTamqcMWYP6qqTsCKvggSgOAJHewc\nr2RZr8j1uGrSSfmIzgdPPQSPQbx1BSy9qQQW2CqY2/xRJ/VjQ7vmaFfWatPEZtCt4B2BSDDznCrs\n+CVE3wbCzsHIy2jodKTm/PYP808C/u7yFi8Ejmjz7Bhs+wkkloGGAT+0/gPt/4MOAqfpi6wp2gcU\nKvVDk63otl/Alsud/Qk2fQFteSDzothbEEsLauAkz7be7/SntSG+IVBzEU5+mg/wOr9Xn4P4RmZe\nHH4C4ktSQQ2c/RHCsOM3NlfUZLAaW4Xr6qofKxeu4YX/vMr2TTs48NhxHDbtIPzB1BppTTekglYc\ndm5913o36h2GBFPL10dfTws86TwQexO8J7U7I1Wnov6JqT41hcDhiM9lVl7kOZypVG1oBBIfgqV7\nmBQLbH1QttSP+x95kGeve4J4LE4yqbw/ezHP/vMlvnvLFfgDYWdzlbZb4WkEWu+FnYFNanBqXYnM\n68QDVGUtk/iGga+zJYuyNTCSqXca47CmaB/Q2FDv2hxNT9iNhGM8cftzRCMxkkmnWRdpjbJ+5UZe\nuv91SO4AyfI9mNy6+/fgZLL+WQXy3OsjdDIQbH9c+oN3VPvjps+ywNZHtU39WP7uSjwC8RGZgSMa\njvHGk2+DdwjutSIPBA7c9Ul8I6D2MsCfmhmQ+ul/LeIJ5Vfo4OTUlKoA4Es9uxb6X2XLFZkM1hTt\nQzpK/QhVB1m/uJq6EdvanavuX42ID625BJr+zO5+Lq8zpanqnIzrJXQ8Gjjc6Y+TAPjH40wnzo+I\nB/p/F40vhdgCZzPl4KR2I7PGWGDrIzpL/Ri133D6D+oHsp34iCC+1c4GxcGqAMedfRSQCljeIdBy\nLyQ3gH88VH3adTqTeKoheES744UgvjHgG9OlezS53VnWSCPOFCzfnj1SNlMaLLD1MdlqbYJwxQ0X\n87uv/IXqwZsI1giJeJyTPnccHz9ydzqI+PeHAfv3ZpHzppE5sOM3qU9J4B/OPqU1X7QmbIWywNaH\ndJb6MXTPwfz8oatZ8u7z1DX5GDNhL/rV1xahpNk5+WrqNEtzuT4Zhqb/R7s0kchTEDgcAuMLXkZT\nfBbY+ri2qR8eEcbufxzh+CL69S+doKbJVmi5BcLPAnHU93Go/a/Om5Sxt3AdI9MIRJ+1wFahbFS0\nj8kl9aPYVCNo5BU0/JyzxwHA9v/eFdQAiL8H265GE1s6e1oHp9wWrDSVwGpsxjVhF5z5o6Orw85U\nJt0K/okQPKogI5zZaPQd2PE/gIAqkEBD0yG+mPbJwTFovgUNHAr+A519SdvyTwBNtD8uwVTOnalE\nFtj6qBWLFzFq6P0QnQviQZPHonohIs7sgJBvLC0t/4LoX9m9zd0cCD+MDrjONbipal6d8aoRJ6i1\nnZIVfgT3HLq4M9IZnQMk0OqzkOqzM64QTxVa+w1ouh6n9pYA/BA8BvwHujzTVAILbH3QmPpqFi/7\nOQxcCyRAoVGeYsmSlTTu81MEQYlC+N8s0RCNVTv35olAfBWEn4Sqabuep+GnoeVOSG5BPfVQ/Vmk\nO/sPROfhvmZpKrC6SrJrwn3Lvah//3YrfUjoSNS/L0RfcvrWAgcjvn26Xj5TNqyPrS+KvggaZsWK\nmrSDcdD1Tt8VQHwxIdfvvQhEXtz1ScNPQ9MMSKb6upKboWkGGn6u6+XSKGjc/ZynDmcFkI5EIDzL\n9Yx4ByFVpyPVZ1tQ6wMssPVBGl9KY906lxNJlq5dkvoQJGvHu1Tv/r3lH7RfcSMCLXd1vWCJNbRb\nInxnWWouhdDU1Ls9ZP3T1dauv9dUHAtsfZF3T3ZOJl+xov+uw42hVvCk5pD6xoD0I4SHJa27r0GC\nUDUVSOWUJbPsEZvc2KUiaWILtN7nftI3GgKHIbUXI4PugPrbce1zkyAEju7Se01lssDWB0nwOJAA\njXVpq3LgA08D+PZhybrNCIL0/4HTBJQgzpJDPghOB/8hznNEdgfCtrzZlyN3FXubrF2+3uEZgxLi\nqYbaL+NMhk/9CUsIfB+D4FFde6+pSDZ40AeJpwoG/BJt/hOwGvBBYBJS+180Sv9dqR/iHQUD/4KE\nH2dpspkxAw9pPy+0+kJoupHM5mgAqi7stBya3AHRV1OjoAIiLq1fyWz67jwaOhH17esMZGgzBCY5\ntbocZySYymaBrY8S3zBkwE+R6CZWbhP2HONewxK8VIWmE44vcp/sHjoGRaD1Lkisd2pqVRchoSM7\nfL9G34Dtv3IGQTVG9lFPP4SmZPln2BNqL+nwPaZvssDWx+0zeFDWVT/GDB0IJJFUf1a2DV8k9AkI\nfSLndzr5ar8GIlnGJwSnDzAONZ9HfHvn/GxjwAKbSclY9UNb0dZ70E1P48zL3IdgzeVE2mb+d1f0\nHTreY9sHVadB1amIp19h3mn6FOuQMDQ2ZG7qotuvg+g8lm5KTYKPL0a3X4smtxVoc+XO5mjGQfwW\n1Ey3WWAzmeIrIL6IxmCbSfEaJtj6YGHe4R/fyQT0IHja9+cZkysLbAZIW/UjuSa1o5Rj6aa0Xdbj\nC9DER3m/SzxVUPs1ss4kkPabJRvTFRbYTCbvnrv2DG0Mbm9zUtHYmwVpjkroKBh4I4ROAapxunv9\n4N0DBvzMCX7GdJMNHpgMK1f4GDlwpLMBccrSTXWMGeQk84akv9uWxd0i3gao/aKzSczOmqB3mC3X\nbfJmNTazy65BhH7fZ+efRmatLYCkMvsLM4jgEBHEt4fzY0HNFIAFNtPOqhUeqLmEXft3Is7voZMR\n31hCvrEdP8CYIrOmqMmwc8MXCZ0C/oPQyAtAhGWR8TQOmlDs4hmTEwtsZSiRSDJv1lu8NvMNfH4v\nR50xifHH7FfQZpyTsDsCqT6P6rUf8cq779E0/F32P3pffD4fId9Ylm5f5DoTwZhis8BWZlSVm75z\nG4vnLiXS6nTjv//6YiadcgjnX3VmQd6xs9aWTCa57Yd389Zz7xId5ecZj4fQz2Jc+MPPMGFyee0t\navqWgvSxichUEXlfRBaLyFWFeKZx9/7rizOCGkCkNcqrD81lrcvuU/l45cE5vP3Ce8QiMWRxC9Fw\nlJbtLfzlyjv4wan/w4alNQUdRDCVrbM4ISKTRWSbiLyZ+vlRrve2lXdgExEvcCMwDRgHnC8i4zq+\ny3TXglc+yAhquynvv764oO965j+vEW3zrnh9FarKlnXbuOHyvxCPFWj+qKloXYgTL6jqQamfn3Xx\n3l0K0RSdBCxW1aWpQtwNnAG8W4Bn9wmJRJInbn2WZ//5Iq1NYUaP34uzvn0aoz42vN21NQOq8QW8\nxKOZW8p5vF6q+xcuqbWxoZ54LPMd/hVRYnsGdpc7nmDRvA/xHe4rSF+basxZ2lv6lVzah8bXQWKF\nk2fnG1Xs4pSjfOJEl+8tRGAbAaxM+7wKOLztRSJyKXApwPA9RhTgtZXj7v+9n9mPziMadjL+F81d\nyv99+SauvvPrDBnVkHHtpGkH89jNT+FsI7ebCEw4rrAV5XFHfYyX1m0l3uS+j0AykWD7mjrXc12h\nGofmW5xFI1Hw1KDVlyChY/J+dkHKtuMGiL4O4gNNoL59oP+1FT87IhKLd2Uj7QYRmZP2eYaqzkj7\nnFOcAI4SkbdxVkD9rqou6MK9u/RaHpuqzlDViao6sb7eJjjvtGNLE6/NnLsrqO0Uj8SYdfvz7a6v\nH1bHxdddQLAmSCj1U1NXwxW/+yLB6mBBy3b2+ZNpGF6PP7D7+8+/Ikq83vkPWhDGHjoGyDNht/kv\nEH4KZyOXOCS3QdONaPStPEpfIK33QXQ2kKpNEoX4B9B8U7FL1uOCHh+NNYNz+gE27vzvO/Uzo7Pn\nu5gH7KmqBwK/B+7vbtkLUWNbDaTXzUemjpkcrF++Eb/fRzyS2VeVSCRZ/u5K13smHDeOXz3xQ5a8\ntRyvz8uYA/fC6yv8d1Qw5Oein57DB3OW8PStz9C8tZlkwlmVI1AV4OAp4xk+ZigwlHB8UbfeoclW\nCD9L+92potByDwSKnDsXfoz2u3DFIfIKWvtV142jjatO44Sqbk/7faaI/FFEGnK5t61CBLbZwFgR\nGZ162XnABQV4bp/QMHIQMZcOeI/Xw4ixe2S9zx/0s9+knt8fc9+hDXgP9zDlM4fz8gOzeX3mPMJ7\neDno+AOYfkYBVuBIbnNW81CXbfeShR3l7Za2u9LvkgRNgAW2XHUaJ0RkGLBOVVVEJuG0KDcBWzu7\nt628A5uqxkXkq8DjOHui/S3VLjY5GNDQjwmTD+CdZxcQjez+j9vn93Hy544rYsky+fw+jj3rSI49\ny9nLYEnzBjxpHfzdTtj1DiLrarqlsLGxfwJEX6PdGua+UYgnVJQilaNscUJELkudvwk4C7hcROJA\nK3CeqirQ5RhTkARdVZ0JzCzEs/qiz//kbO7/fS0v3vs6sUiM4Y1DOe+qT7HHmKHFLtouGUuHpyxZ\nt5nGofVZ7siNiB+tPg+a7wIiaSeCUH1+Xs8uiOrPQ2w+aASnuexzapg1Xyl2ycqOW5xIBbSdv/8B\n+EOu93bEZh6UAJ/fx1nfPo3PfOtUkgktSH/Zu68u4oV/v0K4JcLEkw9i0vSDMwYBumLnTISMYzWD\nd23Tly7bhi8dkarTUKmD8L8hsQV8Y6HmQmcXqg5ocgeEH4XoG+BtgNDpiL+wE/TFNxSt+73T1xZ/\nH7yjIHQK4uvivqmmV1lgKyEigteXf/7Wgzc9wTN3vrArkXfZ2yt4+YHZfPsv/4XX57KDeo7cam3p\nQr6x3R5EkNAx0IX0Dk1uh63fhuQOIOYEncgctPZyJHRst8qQtWzeAVBzbkGfaXqWLVtUYbZu2M6T\ntz+XMTshGo6yZvFa3nhqfref23bDF9i5TV/7PKdemWbVch8kt5M5mhqB5hlOoq/p0yywVZhFc5fi\n9bevlUVao7z9fP6TQVZ2Mh+119Zqi80F1+0AFRKWbdTXWWArI84AUceq+1e5jjF6vB5qB9bk9X63\nWhvQlez0wpEsW/NpHKS2d8tiSo71sZW4RDzBw3+exXP3vEy4JcLIfYdz7vfOoPGgvV2v32/SPviD\nPsLNkYzjPp+XT5zZ4SyUbnEbROiVtdqqToMdS8kYScULvn2cvRRMn2Y1thJ313X38sw/XnQClcKq\n99fw+6/+ldWL3bfB8/q8fP2PX6ZuyIDUtKsQwaoAn7xkCvf8+gG+ftS1fP+TP+exW54hkehs4+L2\ndm3T10Zv19okeARUfxrwg1QDQfCNhn5X9mo5TGmyGlsJ27GlidmPv9luulU8GuOJW5/l4l+c53rf\niH324BcPX82HC1YSbY3Sr76W31zyRyItTu1mx6YmHvvrU2xeu4ULrvk04DRzl761nC3rtrHX/iMZ\nPDL3+byFTP3oCqk+Gw1Ng/gy8NTZqhtmFwtsJWzT6i34fN52gS2ZVNYsXtvhvR6PMGa8kwd253/f\nSzySOVIYDcd49ZG5nHb5ySQTyvWX/Zmt67eDQDKe4JATD+TCH52N1+teqe/J1I+uEE8tBMb3+HtM\nebGmaAkbPHIQiXii3XGP18NIl7Xasln+7krXZqff72P98o389Zq72LhyE5GWCJHmCLFInDeems+L\n973m+rySS/0wpg0LbCWspq6aw085lEAoc6K1P+Dj5C8cn/NzRuwzDI9LzSsWixOqDfHh/BXtAl80\nHOW5f77S4XNLJvXDmDYssJW4c7//KU7+wmRqBtTg9XoYPX5PvnnTpeyRagZGWqPt1nJr66TPTcbv\nz+x18Af9TJh8AKGaYMZk9nTR1ojrcSix1A9j2rA+thLn9XqY/qUTOfasI3n2ny/z/uzFPPPPlzl4\nygHM+vtzLJ+/EhA+dlgjF/7oLOqGDGj3jOGNQ7ni95dw9y/vY+2S9fiDfo4+cxJnfn0aXp+X2oE1\nbP5oa+Z7/V4OmnJAl8patNQPY9qQXJI+C238ARP03ntsMZBcbdu4g+suuJ7WpjDxaNxZ5UfZ/b84\n/W4DhwzgJ/dd2eEk+kQ8gcfrydhT4IM5S/nTN28hHk+QiCcIhALU1lVz1R1fp7au46TeJRs3Zwwi\n7Axs6at+hOMW2HpbfWj/uao6MZ9n7H/gBL3nkSdyuvaAPYfl/b5CshpbGXj05idp2dZKIpEaSNj5\nXZT2nZRMJGna3sKCl9/nwGM/nvVZbpPg9504hmv/+S1e+M9rbFi9kY8d2sik6YcQqun6UuPFSv0w\nJp31sZWBBS+/vzuodSARjbNx9aZuvaNhRD1nfn0ak885mhfufZXvTP4x3578Yx648TES8Y4TeW0Q\nwZQaC2xloLp/dU7Xef2+DpcT78zK99fwx2/cwprF6wCINEd45h8vcvcvs++pYakfphRZYCsDJ3z2\nGIJVgfYn0gYzfQEfw0YPZt/UrlHd8dgtTxNzSeR9feZcmre1dHiv1dpMKbHAVgYOm3oQx55zFP6g\nj6raEIGQnz0/PpLDpx9CTf9qagfWcNzZR/LNP12a10bDaxZ95LqCiNfnZdPaLVnvs9QPU2ps8KAM\niAhnfm0aJ190HCs/WEPdkAEM29u9I761KczLD8xm0dylDN5zEMeedWTO8z5H7TeCDSs3kkxmBrdE\nPEHDCPfgFY/FefLvz/PkM29CNMYhJ45n+pdPpLGfpX6Y4rHAVkZq6qo73HJv++Ym/vfC39GyvYVo\nOIbX5+XFf7/GZdd/gY9NbOz0+VMvmcI7z71LJLx79d1AKMBRnzqM6n7uu57/6Vu3sfiNZcQiMWIN\n1bzw71eZ/+JCfnD3t4DCbPhiTFdZU7SCPPqXJ9mxuWnXTIREPEEkHOXvP/1XTotUDm8cytdvupTR\n4/fE6/PSr76GaV+awlnfPs31+uXvrmLJmx9m9MvFYwm2b9rBvKfe2blDeDs9NYigiQ3o9l+iG89F\nN12INt2MJrPtC2oqmdXYKshbz73rOml+x+YmtqzbRv2wuk6fMfqAUXzvlis6vS6ZVF6bOY9ErM37\nAgEiLVGWvvkhk6Yd0u6+nlr1Q5MtsPV7oDsAdTZgDj8B8WXogF/k1fdoyo/V2CpIsNpl5BTQZNJ9\nVLWbIi0Rfv2F3/PSfa9n5Nf5Nzojp/6gn4ZRziq2vZb6EXk6tfdnes00DvGlkFhc2HeZkmeBrYJM\nPucofMHMSrjX56Xx4NHUDMgtFy4XD8+YxerF69qlhuwkoQBHnHpo1vt7JPUj3naZ8PRzKwv/PlPS\nLLBViHgsznuvLSKZNktARBg8ahAX/9x9pd3umv3oG86cVRcja2u44NrPtJtj2uOpH969ALdcPwFv\n95OWTXmyPrYy5Czj/SFrl65nyF6DGXvIaB6/7VkWvrqIZNq6ah6PMHSvwfQflGVHJxcffbiBR/48\ni2XvLGfQiHqmXTKF/Q7PrGElE+4DEV6vl+/f9jVWtzZnHO+VVT9CU6D136DRtIM+J6j59ivMO0zZ\nsMBWZsLNEW74ygw+WroeTSri9dAwsp6mTU1E2zQNE4kk819aSDQca7dYpZs1S9bx64tvJBqOokll\n80dbWb5gJRf84DNMmnrwrusOPelAXrr/deJpAwciwp7jRhKqCdJYE2SJy9LhPZn6IZ5+6ID/haY/\nQfw9wAvBI6Emv6RlU54ssJWZe294hNWL1hKP7g4qH3U0nSlJaqS088D24I2PEW2NZqSGRMMx/v2b\nhzj0pAk8/DRAAAAPDklEQVS79j847SufZOHsJWxdv41IS4RgVQBf0M/nfnJ21mf3xqof4hsBdb9A\nNQmIBbQ+zAJbmZn92JsZQQ0gEUsgIni8noymKMCwMUOoqg3l9Oyl7yx3zXeLtEbZsbmJusH9Aaju\nV8UP7v4mbz//HisXrqZhRD2HnnQgwerMZY6KteGLiHUd93UW2MpMIu7eaY9AbV014eYI0XAMX8CH\n1+flwh+elfOzBzT0p2lLs8sZparNzAOvz8vBUw7g4Cyr7DY21LNkY+aAgZP6saFdc9TWajOFZl9t\nZWb/o/bD48lsYolH2G/SPvz4P9/ltMs/ycEnjOekz03mx//+Dg3D61n34QbisSwBMc3US6YQCGWO\nLAaCfiZNP5RgDn10bmzVD1MMVmMrM2d/93SWvr2caEuUSDhKIBTAH/Jz3lVnUlVbxQmfPYYTPnsM\nkZYIt/34Hua/tBCvz4vHI3z6m6dw9KcmZX32oScdyNZ123h4xizA6Zs75KQJnPO907tVVrdaG9j8\nUdPzLLCVmfphdfz0vit5/dF5rFy4huH7DOOIUw+hqjazqXjrj+5hwcsLiUfju3LO/vWbhxg4bCDj\njsheSzrhwmM49pwj2bxmC/0bats9N1+24YvpDXkFNhH5NXAaEAWWABer6taO7zL5CtUEOfasI7Oe\n37GlaVdQSxcNR5l16zMdBjZw9i0dmmVZpK5qbKjv9dQPY/LtY5sFHKCqBwIfAFfnXySTrx2bm/H6\n22/aArB5XfG/d3p71Q/T9+QV2FT1CVXdWS14FRiZf5FMvgaPrM+cC57i9XoYm8fS4fmwQQTTmwo5\nKnoJ8GgBn2e6yR/0c8YVnySYNsLp8XoIVAeYdskJvV4e2/DF9LZO+9hE5ElgmMupa1X1gdQ11wJx\n4M4OnnMpcCnA8D1GdKuwJneTzz2aQcPreeLWZ9m6YRv7Tmxk2hdPYMXC1dz+k3to2dHKQVMOYMr5\nRxd8gCCbYiXsmr6n08Cmqid2dF5EvgCcCpygHSzTqqozgBng7ATftWKa7hh/zMcZf8zuzZMf+MOj\nPPvPl4m0OhPF16/YwOsz53HNnd9oN2ug0Cz1w/SmvJqiIjIVuBI4XVU73p/NFNX2TTt46q4XdwU1\ngFgkzrYNO3jlobm9Vo70vja3QQQn9cOaoyY/+fax/QHoB8wSkTdF5KYClMn0gGXzV+JzGSmNhqMs\neGlhr5TBtukzvSWvPDZVzb5lkikp/Qf1w62jwOP1MDCHvRB6Sm+s+mH6Hpsr2kfsvf9I6gb3w+PN\n/L/c5/dy3NnZk317glvqR3qtzVI/TL4ssPURIsLXbvwSw/cZRiDoJ1QTpKpfiIt+ci4jxvbe0tnZ\nUj/cWF+b6S6bK9qH1A8byDV3foP1KzcSbo4wYp9heH3uMxRyFY/F8fq8XV7U0VI/+p7UYOMNgBe4\nWVX/t835zwLfBwTYAVyuqm+lzn2YOpYA4qo6saN3WWDrg4aktsbLx/wX3+Oe3zzEpjWbCVYFOP78\nTzD9yyfuWmW3I7mmfmhyO4vXP0ZjcBn4RkNoOuIdlHfZTe8TES9wI3ASsAqYLSIPquq7aZctA45T\n1S0iMg0nPezwtPPHq+rGXN5nTVHTZYvfWMbN37+Tjas2oUkl3BzhqTte4P4bHunSczpK/dDESoJN\nv4XIsxB7E1ofgq1fR+PLC/GPYHrfJGCxqi5V1ShwN3BG+gWq+rKqbkl9zGuKptXYTJc9MuPJdhvH\nRMNRXvjPa5x6+Sdz2py5s1qbNt8M2gLsXOo8DhqHphlQ998F+KcwnYlF453O8U3TICJz0j7PSCXl\n7zQCSN/gdRWZtbG2vkjmFE0FnhSRBPDnNs9uxwKb6bJ1K9w79cXrYdvG7d1u6makfsR2t1CWtPan\nsWq78yH+HqpqG7X0gqDPlzX30MXGzvq9ciUix+MEtk+kHf6Eqq4WkSE4ebMLVfX5bM+wpqjpshFj\n98A1rqgycMiALj0ra+qHOLW+ULs/0aAFtfK0GhiV9nlk6lgGETkQuBk4Q1U37TyuqqtT/7seuA+n\naZuVBTZDLBJj3YcbaNnRmtP1p/7XSfiDmc3NYCjAiZ87Dn8w970ROkz9CJ5I+paBS1r7O59DU3J+\nvikps4GxIjJaRALAecCD6ReIyJ7AvcBFqvpB2vEaEem383fgZGB+Ry+zpmgf99QdL2TscXDoyRO4\n4JpP4w9k/9PY6+Mj+dqNX+Le6x9m1QdrqK2r4eSLj+fYzxzRrTK4pX5IzYVocg1E3yYkPsIaBX8j\n1Hy+W+8wxaWqcRH5KvA4TrrH31R1gYhcljp/E/AjYBDwx1StfGdax1DgvtQxH3CXqj7W0fsssPVh\nc594i4dueoJoePfE+Hmz3sbn9/HZaz/d4b2NE/bie7dckXcZOt6m71o0sQYSqxDiLJPBjJHOByZM\naVLVmcDMNsduSvv9S8CXXO5bCkzoyrusKdqHPfa3pzOCGjjN0tdnziUajmW5q2dkG30T73AkMImq\nwFG9Wh5T3iyw9WHbNu3Ieq61Kdxr5bBVP0yhWWDrw8aM38t1dDNYHQSBh/70OL+55I/c+sO7Wblw\nTa+WzeaPmnxYH1uFSCSSzH9xIYvmLmHA4P4cPv0Q+g/q1+E9p18xlffnLCEWjpJMOmsaBUJ+pn/5\nRK4777e0NkeIR+Msm7+CN5+Zz8W/uIAJk8f1SPlz2abP5o+aXFlgqwCxSIzrL5vBmiXriLRE8Ad9\nzJzxJF+54RLGHjI6633DG4fy/du/xiMzZrHsneUMGl7PtC+ewLwn36Fle5hEIgGAJpVoOMZd1/2H\nA47ZL6f5oIVga7WZ7rLAVgGe/9errP5g7a5pTrGIsyPi3665i/+eeQ0eT/aE1mF7D+aL112QceyO\nn/1rV1BLF2mJsmXtVhpG9tweBbbqhykE62OrAK89Oq/d3E2AcEuEtUs/6vLzqgdUux5PJhNU1Ya6\n/Lxc2TZ9plAssFUAr8/9/0ZNKl5f1yvlJ3z22HYT2b0+Lx87bCw1de5Br5Bsc2WTLwtsFaCqn/u+\noLUDqxm6V9cnpE+adhDHnnMUvoCPqtoQgaCfvQ8YxRd+dm6+Re2UpX6YQrA+tjLXvL2VxfOWup7b\ne/89uzVhXEQ482vTOPmi41i1aA11QwYwdK/idda7DSI42/QtskEE48pqbGVu3fIN+APuE8+zLS+U\nq5q6aj522D69HtQaG+o73fDFmI5YYCtz9cPqiEXj7Y6LCMOKWMsqNEvYNV1hga3M1Q3uz/5H7Yc/\nmNmr4A/6OOnzk4tTqAKxQQTTXRbYKsDFvziPw6YejC/gw+vzMmiPgXz5Vxex534jil20brPUD5MP\nGzyoAIGQnwt/eBbnfv9TRFujVPevqphVZi1h13SH1dgqiD/go2ZAdcUENUv9MN1lgc2UFbdBBCf1\nw5qjZjcLbKakWeqH6Q4LbKbsWOqH6YwFNlMWLPXDdIUFNlPyLPXDdJUFNlM2rNZmcmWBzZQFS/0w\nXWGBzZQtS/0w2RQksInId0RERaTri38ZkyNL/TC5yjuwicgo4GRgRf7FMaZrLPXDuClEje23wJWA\nFuBZxnTKBhFMZ/IKbCJyBrBaVd/K4dpLRWSOiMzZvHlTPq81fZilfphcdLq6h4g8CQxzOXUtcA1O\nM7RTqjoDmAEw/oAJVrszebFVP0xHOg1sqnqi23ERGQ+MBt5KrSYxEpgnIpNUtet7vhmTo8aGepZs\nbF9DS9813vRt3W6Kquo7qjpEVfdW1b2BVcAhFtRMMVjqh0lneWymLFnqh+lIwQJbqua2sVDPM6ar\nLPXD7GQ1NlPWLPXDuLHAZsqWpX6YbCywmbJntTbTlgU2U9Zs1Q/jxgKbqTiW+mEssJmyZ6kfpi0L\nbKYiWepH32aBzVQMG0QwO1lgMxXBUj9MOgtspqJYrc2ABTZTQbqS+mG1tspmgc1UtGypH6ayWWAz\nFcVSPwxYYDN9gKV+lAYRmSoi74vIYhG5yuW8iMjvUuffFpFDcr23LQtspiJ1Vmuz5mjvEhEvcCMw\nDRgHnC8i49pcNg0Ym/q5FPhTF+7NYIHNVJxsqR9urNbWayYBi1V1qapGgbuBM9pccwZwuzpeBepE\nZI8c783Q6Z4HPWH+grc37rv/yOU98OgGoFwWuyynskJ5lbecygo9U9698n3A/AVvP77v/iNz3QQ9\nJCJz0j7PSG3gtNMIYGXa51XA4W2e4XbNiBzvzVCUwKaq7l+feRKROao6sSeeXWjlVFYor/KWU1mh\ndMurqlOLXYbuKkpgM8b0OauBUWmfR6aO5XKNP4d7M1gfmzGmN8wGxorIaBEJAOcBD7a55kHgc6nR\n0SOAbaq6Nsd7M1RajW1G55eUjHIqK5RXecuprFB+5e0yVY2LyFeBxwEv8DdVXSAil6XO3wTMBKYD\ni4EW4OKO7u3ofaJqm7IbYyqLNUWNMRXHApsxpuJUZGATke+IiIpIrjk4RSEivxaRhanpI/eJSF2x\ny9RWV6eyFJOIjBKRZ0TkXRFZICLfKHaZOiMiXhF5Q0QeLnZZKknFBTYRGQWcDKwodllyMAs4QFUP\nBD4Ari5yeTJ0ZypLkcWB76jqOOAI4IoSLy/AN4D3il2ISlNxgQ34LXAlUPKjIqr6hKrGUx9fxcnP\nKSVdnspSTKq6VlXnpX7fgRMwRhS3VNmJyEjgFODmYpel0lRUYBORM4DVqvpWscvSDZcAjxa7EG1k\nm+JS8kRkb+Bg4LXilqRD1+N8CSeLXZBKU3Z5bCLyJDDM5dS1wDU4zdCS0VF5VfWB1DXX4jSj7uzN\nslUqEakF/gN8U1W3F7s8bkTkVGC9qs4VkcnFLk+lKbvApqonuh0XkfHAaOAtEQGnWTdPRCap6ke9\nWMQM2cq7k4h8ATgVOEFLL6kwl2kwJUVE/DhB7U5VvbfY5enA0cDpIjIdCAH9ReQOVb2wyOWqCBWb\noCsiHwITVbVkV3kQkanA/wHHqWrJrZ8jIj6cQY0TcALabOCCzrK+i0Wcb7TbgM2q+s1ilydXqRrb\nd1X11GKXpVJUVB9bGfoD0A+YJSJvishNxS5QutTAxs6pLO8B95RqUEs5GrgImJL69/lmqkZk+piK\nrbEZY/ouq7EZYyqOBTZjTMWxwGaMqTgW2IwxFccCmzGm4lhgM8ZUHAtsxpiK8/8B86wDQrw+WR0A\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4476ff41d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_train, y_train = create_toy_data()\n",
    "x1_test, x2_test = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))\n",
    "x_test = np.array([x1_test, x2_test]).reshape(2, -1).T\n",
    "\n",
    "feature = PolynomialFeatures(degree=1)\n",
    "X_train = feature.transform(x_train)\n",
    "X_test = feature.transform(x_test)\n",
    "\n",
    "model = LogisticRegressor()\n",
    "model.fit(X_train, y_train)\n",
    "y = model.proba(X_test)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)\n",
    "plt.contourf(x1_test, x2_test, y.reshape(100, 100), np.linspace(0, 1, 5), alpha=0.2)\n",
    "plt.colorbar()\n",
    "plt.xlim(-5, 5)\n",
    "plt.ylim(-5, 5)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3.4 softmax回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们对于多分类的⽣成式模型的讨论中，后验概率由特征变量的线性函数的softmax变换给出，即：\n",
    "$$p(C_k|\\phi)=y_k{\\phi}=\\frac{exp(a_k)}{\\sum_jexp(a_k)}$$\n",
    "其中，$a_k=w^T_k\\phi$\n",
    "$$E(w_1,...,w_k)=-lnp(T|w_1,...,w_k)=-\\sum_{n=1}^N\\sum_{k=1}^Kt_{nk}lny_{nk}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD8CAYAAACLgjpEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl8XFXd/9/f2dOsTdJ9Lw1dKd1oCxQoUAoUEFlEVBBQ\nxMoPH/QRFUVRH0FFfQS0KuCDgoqCAkqFsrRAS1ta6EJLW7okpXvTJU2zZ5Y79/z+uJN0JjOTTJKZ\nzCQ579drXpm5y7lnZnI/c873fBdRSqHRaDSxsKW7AxqNJnPRAqHRaOKiBUKj0cRFC4RGo4mLFgiN\nRhMXLRAajSYuCQuEiPxRRI6JyNawbYUislRESkN/+6ammxqNJh20ZwTxFHBZi233Am8qpUqAN0Ov\nNRpND0Ha4yglIiOBl5VSk0KvdwJzlVLlIjIIWK6UGpuKjmo0mq7H0cnzByilykPPjwAD4h0oIncA\ndwD0yc6afvrY0Z28tEajicemjdsqlFL9OttOZwWiGaWUEpG4wxGl1BPAEwBTp09Sb737j2RdWqPR\ntKDQM3FfMtrp7CrG0dDUgtDfY53vkkajyRQ6KxCLgVtCz28BXupkexqNJoNozzLn34E1wFgROSgi\nXwR+BlwiIqXAvNBrjUbTQ0jYBqGU+kycXRcnqS8ajSbD0J6UGo0mLlogNBpNXLRAaDSauGiB0Gg0\ncdECodFo4qIFQqPRxEULhEajiYsWCI1GExctEBqNJi5aIDQaTVy0QGg0mrhogdBoNHHRAqHRaOKi\nBUKj0cRFC4RGo4mLFgiNRhMXLRAajSYuWiA0Gk1ctEBoNJq4aIHQaDRxSYpAiMjXRWSbiGwVkb+L\niCcZ7Wo0mvTSaYEQkSHAfwEzQjU77cCNnW1Xo9Gkn2RNMRxAlog4gD7A4SS1q9Fo0kinBUIpdQj4\nJbAfKAeqlVJvtDxORO4QkfUisr7ieGVnL6vRaLqAZEwx+gJXA6OAwUC2iNzU8jil1BNKqRlKqRnF\n/Qo7e1mNRtMFJGOKMQ/Yo5Q6rpQKAC8C5yShXY1Gk2aSIRD7gdki0kdEBKsU3/YktKvRaNJMMmwQ\n7wHPAxuBLaE2n+hsuxqNJv0kXLy3NZRSPwB+kIy2NBpN5qA9KTUaTVy0QGg0mrhogdBoNHHRAqHR\naOKiBUKj0cRFC4RGo4mLFgiNRhMXLRAajSYuWiA0Gk1ctEBoNJq4aIHQaDRx0QKh0WjiogVCo9HE\nRQuERqOJixYIjUYTFy0QGo0mLlogNBpNXLRAaDSauGiB0Gg0cUmLQHgNIx2X1Wg07SRZxXsLROR5\nEdkhIttF5Oy2zimr1NW1NJpMJylZrYFHgdeUUteLiAurPmdcbOIGIkViTKGutqXRZBqdFggRyQfO\nB24FUEr5AX9b53nsJRGvyypLtUhoNBlGMkYQo4DjwJ9E5ExgA3C3Uqo+/CARuQO4A2DQ0KFRjXjs\nJZRVlkZs04Kh0aSXZNggHMA04PdKqalAPXBvy4PCi/f2LYp943vsJc0P0HYKjSbdJEMgDgIHQyX4\nwCrDN62zjWqR0GjSTzJqcx4BDojI2NCmi4GPOtsuaJHQaNJNslYxvgo8E1rB+Bi4LUnt4rGX4A2W\n6hUPjSYNJKt47yZgRjLaikX4ikeTWGiR0GhST7dztdbTDo2m6+h2AgFaJDSarqJbCgRokdBouoJk\nGSnTgjZgajSppduOIJrQzlUaTero9gIRjhYJjSa59CiBAC0SGk0y6dY2iHho20TiKKUgsBmC5eAY\nCY5xiEi6u6XJEHqkQIB2rkoEFayGmvvAPAnKALGDfRgq74eILavr+2McAPM42Ecidv1dZQI9ViDC\nCR9RaJEIo/53EDwKBK3XKgDGHmh4BnJu77JuKLMOah60ri12UAGU5yLI/nLM0YxSAfCtAmMH2AaB\n5yLEltdl/e1N9AqBAC0SLVHKAP9GmsWhGQN873SJQKhAKTQ8C4GtQCC0MbTT+wZ430DZCiHrBvBc\ngoigzFqo/haY1aC8gAsa/4nKfwBxjEp5n3sbPc5I2RragBmOAsw4+1qKRgqu7t8C1d+HwAc0i0Ms\nzEqo/xN4X7ZeNzwLwRMhcQDwg2qEukdT3eVeSa8SCIgUiaZHb0TECY5xMfbYwXVW6jvQ8CQJZCYM\n4YOG51DKBP8aIEZWdOMwyqxJYgc10AsFAiKdq6xUd71TJMj5fyC5EEoijHjAVgB9bkn9tY397Tte\n+UHV0/qs2N6ZHmli0GtsEK3RW/NhimMwqu/vwbcCggfAMQbc5yJNgpHSi+eCqm3H8S6QPuCZBw0v\nEjn6sFnLs7bsZPey16MFIkRvXRYVWx/IurzrL5x1NTT8E/CFbbSB5IWEI9wO4oas6xCxo7KugcAO\nCGwHlLXqYcuF3Lu7tPu9hbQIhN+X2ZW19IpHF5B1DZh14F0SWtoMgmcBZN8M/nXQ8BQEj4AtH7Ku\nA88VQMh2kn8/ytgNxm6w9QPnmYj0ytlyyhGlVNtHJZmS8RPVo396lmGj+3f5tduDN2hNO7RIpA5l\nesE8AbYixOaJ3KeU9ursIIWeiRuUUp3O8pYW2fXYrYHLgY+PpePyCaNXPFKP2DyIY0iUOABaHDKA\ntI3LSvKsX+XuIBI6nFzTW0maQIiIXUQ+EJGXEz2nu4hEE1okNL2NZI4g7ga2t/ckLRIaTeaSFIEQ\nkaHAFcD/deR8LRIaTWaSrGXOR4BvAbnxDggv3jt40JCo/SV5hZTWVEaIRCavcuicE5reQKdHECJy\nJXBMKbWhtePCi/cWFhbFPKYkr7D5AZk/ouipBkwVPI7yb0IFj6e7K5o0k4wRxLnAJ0RkAeAB8kTk\nr0qpmzrTaPiIIpNHEk30BOcqpQJQ+wj414M4rbwMrumQ+3XLQSnp11Pgfx+8y4AAuOeC+zxEdExF\nppCM4r3fUUoNVUqNBG4E3uqsODTRXUYSTXT7kUT93y1xIACqwfrr3wD1z6Toeo9B3SMQ2ACBD6Hu\ncaj5Celw3tPEJuP9U7uzSHQ75yrf60TnZgiAbynK2IOq/jGq8jZU1bdR/g0osx7lXYnyrrASubQD\nZRwA7wpQ4bEYPjC2W2KhyQiSGouhlFoOLE9mm9A9pxtNdKtpR8TNGr69Eaq+Q3MEpVkNNT8DJDQV\nUUAQlbMQ8VyY2LUCWziVPir8Wl7wfwCuM9vff03SyfgRRBPdbSTRRLeadjjGxt4ufYhO7hIEDEs8\n8AIBqHsMZST4/UgOSKzfJwfo/JIZQ7cRCNAikXJybgfJ4tTA0mElkUnYJqDAvzqxQ+NmrbKB+4IE\nr6dJNd1KIKBniESm2ibEMQoKHgHPpeAYD575kP8IOIoTbCEYf5rS8lq2LMj+IpFZoByQ+w3EHnsZ\nXNP1dDuBgO4tEpme6k7s/ZCcLyIFDyI5tyOO/uC54VRaulZxgntmQtdRxgGo+wORiWFsENjWkW5r\nUkS3FAjoviIRTqaKREvEcy70uTlki3ABTnCcFnreFJLtBs/FiGN0Yo02Pk+0XcMP3ldRZkOSeq7p\nLN065Vx3W92IRXfJhylZC1Ce+dYKhi0XEZdV18K3AjDBPcealsRAmfVWrY3gwea8lxi7ibmKIQ4w\nj4JN17jIBLq1QED3FwkVrMDduNyqEmUfgtc9jbLKrhEJpUzrumY9OMchtrihNACIOCBkH1C+NdD4\nHAQrrdGEeGJXwTIOQPV3rdJ++CyjZ8OzYB8KwcMxOuUDW6I2D02q6fYCAd1XJJRxGKq/GTLsBcHY\ni8f3Ht6cmymrHJtSkVDGAaj5UWiZEsutOvtzSNbVbZ/b+CrUP03zFCGwGap2oAp+El3dqm5RKF19\n08leq8SfYyTW9CTGKMLYA67JHXhXmmTTbW0QLemWNomGv4QqRDUZ6kzAj6dhMUqplNknlFJQ8z9W\n1SrVGBIJA+r/jgq0ntJDqSA0/I1o+4EP6v8WeazpDU0lWhIE/1Zi17EwofGlxN+MJqX0iBFEE90t\nZNyqSRmj/J15Eo9tID51NDXh5MaOUKxFS/zgfQ2csW0JTX1DxSmVFyyLfC12ThkxW2ILeWHGyHCu\nMt9w21voMSOIJjI1ZLzW52Xlvr2sObgfnxG6KWw5cY4WEE/qwsmbphWxaCumojUvR1ukGIs4wTmZ\n6H8zJ3jmErs2qAOc01rvg6bL6HECEU6miMTS3WXc9tKLLFq3lkfXvsvn//08Hxw5DJ6rgZb+BS5w\nn4OIK2JrUkXCMS72LzducJ/T6qkiLqu6Fa4We9zQ54boE3LuAnv/kIem0zJSOkZD9k2hpdPw9++w\nRNPzifa9H03K6NECAekXiUM11TyxYR0B08RrGDQaBl7D4KcrV9Bou9DyVsR1ysfANRmyF8ZsK1ki\nIbY+kH0rEX4M4gbHsMTcnLNvO9VvnFYZvZw7ENf06GvZ+0LBIsi9B3JuhbzvQ/5PEHEjWQsg915w\nTgHHcMi6EvIfRuz5nXp/muTRo2wQ8UjnKsfyvXsIquihtIjwfvkh5o78AqrP9RA8BLZ+iD1yic8I\nBllXfogjdbWMLOjLlAFj8JllnbZNSNblKMcY8L1mTStcs0PJWtpODCNih5wvoLJvtpZIbXmtVrYS\nsYFrKhAykHqXoXyvgOm1rpv79TaXWDXpoVcIBKRPJLyGQTBGsJNpmmw+coSPjh9jcG4eF40cTZ4z\nsnhMRX093172OnUBP0YwiMNuZ3BOLj+5eD5ZTutG7kw4uThLwFnS9oHxzhcn2Avad1L9b8G7iuZV\nEO9/wL8aVfBrKz5Dk1H0+ClGOOmYbswcMhS3PVqH/abJ6gP7eH13Gc9s2cwdL/+bPScjpw6/WbeG\nysYGvIaBoRRew2B/dTV//XBT8zFN047SEydYuW8v33trKd9e+jqvlu4iEIxlBEwfKrAHvO8QuUSq\nrNJ7Dc+mq1uaVuhVAgFdLxKT+g9g1pCheBynRMKGNfP3BS3/B38wSKNh8PDad5uPCRgGW44ejbLz\nG8pkxb49Eds89hL+sW0Lj763hg+PHmV7xXH+tGkD97+9DNNMf/o2pfyo2l9B9beAOIWbva/rVHMZ\nSK8TCOhakRAR/vvsc7n33PO5aOQoLj1tDH1crlj+gxyqqaHG6wUs/8J4N0zL7Qeqq3h5px9/8FRk\npNcIsvtkJesPH0zWW+k4db8F33tERm62xACdRTvj6JUCAV0vElMHDebu2edy51mzI0YTkSjsNmtV\nweVwML64f5SbkV2Ec4aNiNi27fgxEDhUE5lHwWsY/POjrUl6Fy16qsyEfvGV2QC+NUTnumxJyHFK\nk1H0WoGA9C2BXjJ6DE5bpJuxDRhX3I9s1ym/gK/OOptct7tZUDwOB/2ys/n85CkR5+a63NhCgVKH\naoqaHwBlJ06wat++pPVdGXtRVd+BE9fDiRtRdb+3XKrjYdYQ26W6BY5R1pKoJqPo9CqGiAwD/gwM\nwBoZP6GUerSz7XYV6VjduG7cBD46fowdFRUopbDZhFyXm6/PPjfiuEG5ufzhqmtYuW8vR+pqGVXQ\nl1lDh+Gw2VBKNUdPzhw8BEeMSMpDNUUMyTvB37Z+yJwRI6L2txdllEP1d8KyRgXAuxyCRyD/Ry2O\nPQbBfWArslyuYw42LI9RJNvyk9BkHNJZw5CIDAIGKaU2ikgusAH4pFLqo3jnnDHpTPXiP5Z06rrJ\nprTGWkHoyiXQ0hMVlFWeoH92DlMHDsJma31At73iGI+vX8eeqpN4HA4uH3M6N51xJg67nXWHDvLA\nyuUxzxtbXMuP5l5MSVHHUrkppaxVhsYXiO0e7YKCXyCOYVYwV92vwbfWyu2ggmArtALDmlcvxHLM\n8lwHztHgPLNVPwpN+yn0TNyglJrR2XY6PYJQSpUD5aHntSKyHRgCxBWITCQdgV4lRcWUFCWW++BA\nVRU/ePstfEFrFcBrGLxSupOT3ka+Pvtcpg8aTI7TSV0geq5fWpHLI++9y50zZuEOTVdi+U0o0wu+\npeBbD/Z88CxAnOOsAK7GxcQWB6wRQvCQ5YnZuDhkkAycCuoyK8BRcuq5cwJk3YA4BiX03jXpI6mO\nUiIyEpgKvBdjX6vFezOBJptEE6VdOO04Xl/H8YZ6huXlk+v2RO1/fsc2AsHIJUJ/MMiq/fu4bcp0\nCjweFs6Yxa/fXxOxmgHWbb3+oIu/uar5ylmzASirLI0QCWU2WrkpghWA31qN9L2Pyr4NGl8EWklG\nqwywD7eee18lOhTcAKMUiv4aFWOiyWySJhAikgO8AHxNKVXTcr9S6gngCbCmGMm6biopyStMuUh4\nDYNfrH6HzUeP4LTZCZhBrigZy61TpkVkaNpbdTL24N5m40hdLQUeD+eNGElRVhbfeWtp1HGGMlm+\nb2+zQLRMdae8yxnjColDM36of4rWlyed4JqMOAY3NRTnOGWNKLRAdCuSMvETy4H/BeAZpdSLyWgz\nUyjJK0zpKsfv3l/LpqNHCJgmDUaAgGnyatkulu6OzK1wWkFhzC8rYJoMyjkVxzC+X/TSaBOGGSkx\nEeHkgW2U1fSJPklsYG9lKuBZALnfOvXaNYWYOSDsAxFbdvx2NBlJpwVCrJ+5J4HtSqlfdb5LmUeq\nRMJnGKw+uD/qxvUFg/xrR6QJ57oJE3G18J9w2e3MHTmafM+pKYmIMHnAwKhbVIAZgwbH7YvHNgCA\nspqWsRUmZF1DdHi3C/K+h+TcEhng1edmq2oWTdtCxXeyv4LyrUPVPY1qfKXdtTw16SEZI4hzgZuB\ni0RkU+ixIAntZhSpEIkGI77zUK0/cs4/JC+fBy6cx7jiYuwi5LpcXD9+InfOiK5D8ZUZs8hxuXHb\nLf8Dj91BvtvD7dNiG7WVUhwPzsWlLKE5JRJiJZB1XwAFPwXXTLD1s8Kz83+IuKzELl7DwAjFfYi9\nH/T9DfS5zjrOcznkPwQNT0Pdw+B9Cer/AicXWlmxNRlNMlYxVhE/r1iPItk2iQK3hwK3h4rGyPRv\nAkweMDD6+kXFPDTvsjbbLczK4q6Zs9ly9CheI8CYwkLmjhzdHAEazq6KCh5a/Q51AT8Lho7jxtHb\nMOw2ymr7M6bADrnfs2whjlGQd2/EuTsqjrPo/bUcqqnBZrMxZ9gIFs44iyxnXkTyGNXwIhh7OeVN\n6bf8Imp/ier7WMxs2JrMoNeEeyeLZIqEiHDnWbN4aPU7+INBFOAQwW138JmJk1m9fy9H6uoY3beQ\nKQMHJXQjrdq3l9+8vxYR6x502x1cNqYkpjjU+rx8f/kyvKEUeC/uLeG1gyOYXFjDXWdPoswcjtQI\ncCqcXAVPgqriSEMu97/9ZvO5wWCQVQf2cdLbwP9cOC/yQr4VxHS1NmsgWA6O+FMfTXrRAtEOvPU+\n1vxnHTveK8Mcnc8Vt13AwJGdE4rpg4fws4vn8+KOjzhcU8P4fgM4b8Rw7l++jAbDIGAYOB0OBufk\n8s1z5rCr8gRZDifTBg7C2cImUV5bw6Pvr+G03CPcVrKFkbnVVPs9LN65g5Ez7406fuX+fQRb2D8a\nDBebKgey/shoLhxlVcnyBksprTjMGNefwL8ZxEmRGWD+4Aks3j+m+dxAMMj2iuMcqqlhSF5Y7sq4\nTlCKXu7tn/FogUiQ+uoGfnbTr6k9WYffG8D+ro0P39nGl358I5PmjOtU26MLi7jnnPOaX9/35htU\neb3Ny5pBw2Bv1UnuXPIf3A4HAgjCD+ZexLjifs3nLft4N6NzK/jB1NV47NbSZD9PA587bROHjz/F\niEG3R1x3X9VJAmb04qnPMKjynkps67GX0Fj/Xcpq9jEmzwBl4BT43JhtlDdks67i1CqHQ2wcqauL\nFAj3JWD8mSj/CHs/xBE9ldJkDlq+E+SNp1dQXVGD32sNlYNBE/ueWp5+8EXMYGt+Au3Daxhsrzge\n5fNgYv3eNuW1bDAC/HjFWwTCrl3r93HDqK247JH98diDDHYuRbVIV7+nqipmHxRECI8yq/H4ywAj\nYpXDYw9yzchdEecGTJPh+S1ySnrmg+sMrAS9DiuBreRGLo9qMhI9gkiQzcu3YgSihcCxt5bNa8qY\nOmdsUq7TntiYoFJsOXqEaYMtz9QZg4YwylYTU/VtYNW0sJ+aEh2sqY7b9tDcsBGAWQtix6NseMWM\nEIkiz6mRhtth5+yhw+mXHenvIOJA5X4XjDKrJof0BfeshPJfatKLHkEkiCe7ZXp6C9M0GV9YnLQl\n0Cynk5KiooSXhXxhI4gZg4dQ7Y+dm9KGgC3ylz3XFfs9OWy2yJwV9gE0/at4lK35AcLGynFkORwU\nZmXxqQmTuHtW7LT5IoI4S5CsqxDPHC0O3QQtEAly4afn4PJEOgvZ7DaGjh1CQf/8Zj+J8EdH+dqs\ncyLyQDjjRHkapskZ/Qec6o/NxrCBX8ag5c3nsgKvJFIQPjlufLOvRBNOm40Lho+MMGiKOK1U9xHO\nUnY85HD6kBt58OL5PPXJ67hh4hnNCW80PQM9xUiQmVdMZffmvaxZvB6FQhAKBubzpYduaj4mPNir\nMzkmBuXm8YerrmH1/r0cratjZEFf3vi4jO0Vx/EaBgI47Xa+OHU6Oe7Im97hmYySb0LDk1YKN/FA\n1icg61MABE0TmwgiwmVjTudwbS1LynbhstkImCbTBg3myzGcr8QzD2XrZwVumRXgPAOyrifLXtyp\nzNqazKbT+SA6Qibmg2iLxrpGfnbTb6g6Xk3AZyACDpeDm75/PWddNjXmOcnMMWGaJhsOH+bdg/vI\ndrmZN/o0Rha0noHJMko6EBE+rjzBb9e9T9nJEzhsNi4cMZrbp8/A43BQ6/NyqLaWfn36UNSnY/ES\n3qDlFalFIjNIVj4ILRAJsuT/lvH6n94m4IsMufZku3lo6f04XbEHY22JRJXXy9G6Wgbl5JLniQ7z\nTgbH6+u469WXm52awJpKTOjXP9qpqRM0iUQTWizSR8YkjOktbHp7W5Q4NHG4tJwRE4fF3BcvpZ0R\nDLJo3VpW7t/XPLy/aNRoFk6f2WZmqfbySulOjBZLsQHTZNuxo+yoOMa44uS4jjfV6IDOFfTRZA7a\nSJkg2XkxQqEBM2jiyW29IlSs5LjPbNnM6gNWJGeDYRAwTd7euyciC7UKFcvp7ChvT9VJjBhtGEpx\n35tL+cvmD5JekyLpFck1aUELRILMvfFc3FnRqxj9R/RjwPC208aFi4RSiiVlu6IyP/mDQf6zaycA\nK/fu5QsvvcBnnn+Wz77wHP/ctqXDN/HphcXxV0KU4j+7drB6f/IyXzehRaL7owUiQc68YAIXfe48\nnG4Hnmw37iwX/YYWsfB/b0m4jSaR2P/xsQh7QDgNAT/rDh3k1+vWUBlyt24I1bd4btuHHer7gpKx\nuOzxU8/7gkFe2rmjQ223RbhIND003QdtpGwntSfr2LftALl9cxg+YWiHQpVLayr55bur2OtoiNo3\nrrgYnxFkT9XJqH0eh4NnrvkUjrCbven7a6sfh2qq+c37a9heURFz/8CcHB6/8pPteRsdQq92dA3J\nMlLqEUQ7ye2bw6Q54xkxcViH8xiU5BVy7fiJZNeoZo9JG1Zo9u1TZ3Csvi7meYZp0hDKWl1eW8sP\nli/jmuee4fp//p2H166m3h8/seyQvHwevGg+Oa7onJAOEWYOHtqh99JeYo0o9Kgic9GrGGni0hGn\nkeVw8Oaejzni9jG6b1+unzCJIXn5DMvPZ0eMX3qP3U6Oy0W938c3l75Gnd+HwhKOVfv2sb+qil9d\nuiCucNltNv7fjJk8/N4aAqH8E06bjRyXi+smTEqo3ye9jRyqrmZATg79snM69N7DVztAr3hkMlog\nupiA32DftgM43U7mjBvOoFwr4Wz4EujnJ0/lhyveijBiuu12PnfGmdhsNt7aswdf0IgoVmUok8N1\ntWyvOMaEfgOIxznDRzIgJ5fFu3ZwrK6OKQMHsaDk9Jip9sMxTZPfrnuP5fv24Apl354+aDD3nD0n\nKs9Ee/HYS7RIZChaILqQD97ayl9+9A8AlII+uR7ufPQ2Gvq7I/wkJvYfwP3nX8hTmz5gf00VRVlZ\n3DhpMnNHWglc9lZVRq2AWG0qDlbXtCoQAKcVFkWV+WuLF7d/xDv792KYZnOS3Q3lh3ly0wYWzpjV\nrrZioUUiM9EC0UEaaht5adFrrH99E0oppl58Btf+1xVkF8T2lzi6v4Knv/8sft+pnAy+Bh+PLvwD\nP3n1u+xprIkQiTMGDOR/L708Zluj+xbi3r8vIpITLEPlsIL8mOd0lv+U7ogSpYBpsuzjj/ny9JlJ\nySupRSLzSFZdjMtEZKeIlInIvW2f0b0JBk0e/tJjrFm8jsY6L956H+8v+YCf37qIoBE7ecyal9YR\nDEZnbwoEDLavLU240nh5bS0byw9HiYMDyHO7eeqDjdzzxqss2bUjIplMZBs1vLN3D9uPH0vYt6Ih\nRkk/AMMMRnlpdga9LJpZJKO6tx34LXAJcBBYJyKLWyve293ZsbaUisMnIxLIBI0gNSfr2LT8I6bP\nOyPqnJrKutjiYVrp7KDtSuMnvY3c88YS6lvcrDagKDubKq+XY/X1AOyvrmL1gf08cNElzb/upmny\nyHvv8u6B/dhDjlPFWX144KJ59M2KPfJpYkK/fmw6ciRq+7D8/E7bIFqiXbYzh2SMIGYCZUqpj5VS\nfuBZ4OoktJuxHCwtx4gRl+Gr93Fo1+GY50w8d1yUJyaAaQYpmTa6+XX4SKJlXoklu3biDa0+hGMT\nobLRGzGq8AWDlJ2s5IPyU/1ZUraLtQcPEDBNvIaB1zA4XFvDL99d1eZ7/uIUK/LTHhIba1nWzldi\nhIYnE+2NmV6SIRBDgANhrw+GtkUgIneIyHoRWV9ZeSIJl00f/YcX4/BE/2q6+7gZMLxfjDNgytyJ\nDB4zMCLpjDvLxXnXzqZocGTYdkleYfMDTk07dp2oiKrCBZbtwVTR272GwYfHTv3qLyndGTU1MYEd\nJyqo8carqWkxvKCA31x2JZeeVsKYwiIuGjWa/51/eZsG0WSgRSJ9dJmRsjsW743H5PPHk52bRcAb\nwGyqKCXgcNiZcnFsfwK7w8bXHv8yaxavY/3rm3FnuZhz7WwmXzC+1WtFTDvyC9h67GhU4JVSCofN\nRrDFze+IipVIAAAWXUlEQVS02Sj0nJo6+OLYRwTBF4zt+h1O/5ycmMlkuoJwA2YTetqRepIxgjgE\nhMc6Dw1t67HYHXbu+eOdjD3rtFMbRQj4Azyy8HEa62L/GjtdDs6//mz++w8Luf3nN5NbmM2x/bFd\nn8NpGknMcBVHuFmDlWa+pKgIR4xgLLvNxvkjRzW/nj10GI4Yqw0FHg/FrSSKafD7CcSJHelKIooN\no0cUXUEyRhDrgBIRGYUlDDcCn01CuxlNQf98Cgf2xeGyY/iDKFPh9wY4VFrOC4+8zE3fuz7uuav+\n9R4v/OplxCaYQZOBI/uz8OFbKeiXF/ecJpG4a/Q0nj+2i/1VVVa5u+EjWDh9JuW1NTy4agV1fqv2\nhNvu4FvnzKEgLAnNpyeewdpDB6j1+fAFgzhEsNtsfH32OTGXKbdXHGPRe2spr6tFEGYPHcadZ80i\nO4a7dlejl0S7hqQEa4WK9T4C2IE/KqUebO347hysFc7Xzv1ehF9DEy63k0dWPxDznN2b9vCbu55s\nrq8BVtj44DED+e4zd7d5zaYMVf2HF+Kw2ZpXI8CaauyrrsI0TUYW9GX3yUpeKd1JZWMjZw0ewrzR\nVhWst/fs5sNjRxmSm8dlp5XQLyfaZfpIXS13v/ZKRNSpw2bj9MIifjrv0jb72VXo4K/YZFRGKaXU\nEqD73/EJYAQMGmq9ZOf3IWhEGwaBuL4QAG/+bVWEOICVdObYvgoO7z7K4NNaN/o12SSO7a+MWgoV\nkeY8lW/v2c3v17/fXPNze8VxXindya/mL2DB6eNYcHrr1cD+s3NnlB+FYZqUnaxkf1UVwwsK4pzZ\ntWjbRGrRnpQJEjRM/r3oVVY+vwbTVLiz3AwYWUz5nmMo89QoTGzCuJklcdupqaiNud3usFF3sj6h\nvrTlL+E3DB7bsC5ixcIfDFLR0MiSsl1cn0Bg1sGaKoIxRpcOEY7W12WMQID2m0gl3S7c2+8NsPQv\nK3jws4/w81sX8e7i9Zhm6hdF/r1oCSufX4PfG8DwG9RX13N0bwWuLBcuj1WHwuVx0ie3Dzd8O74b\nyKQ543C6o3XZCAQZNi7xKteteV5+XHUyZuGdgBlkzYH9CbU/rl+/mFmoAqbJyPzMEYeWaANmculW\nI4igYfLwHY9xePdRAqG5f/nuo+x8v5TbHvhMyq4b8Bus/OfaKHtDMBjEZti46iuXcqisnGFjBzP7\nyun0CeWobKxrZPW/11G68WP6Dy/mguvP4fxPnc2qf71HTWVds7OV2+PiyoXzycppX1breCOJbKcT\nM45tqWUdjXhcPmYsL+/aieH3Nztmuex2zh46LKbNIpPQBszk0a0E4sMV2ziy93izOAD4Gv1sensb\nh8qOMGRMaipFN9Q0xr3hAr4AfQfkM++m8yK2V1fU8rObf01jTSN+XwC7w86qF97jzke/wHee+RrL\nn1vNlhUfkVuUy0WfmcP42WHD5Hofh8qOkF+US/HQ1v/BY4nEsPwCBuTkcqC6KsLr0m23c9XpidUQ\nLfB4eHj+Ap7evJEPjpTjcThYMGYs14yfkND56UaLRHLoVgKxc10ZvoborEkC7N60N2UCkdM3G7vN\nRjxPgB3vlzH9kskR2155Yhl1lfXNzktBI0jQCPLkd57hp6/dxxVfmscVX4quSbH0L+/wymNvYHPY\nCBomw8cPYeEvb4kbJQqxReL+8+Zy//I3qfQ2YgP8psl14ycyox2Zo/rn5PDNc89P+PhMQxswO0+3\nskEU9C/A4YpOvmpz2Mgvzk3Zde12G7OunB57n8NO3wHRIdZbVmyL8mwEqDlRy4p/ronZ1paV23nl\n8aX4fQG89T4CvgB7tx7gD/f+tc0+trRJ9MvJ4XdXfIIHLpzHN86Zw1OfuJZPT5rcWhM9Eu1c1Tm6\nlUDMvmo69haehCKWh+LEcxMbOneUT91zFZ7saBuBw2nn7Kuil5tdfeI7E73y+NKYYdZvPrMSv9cf\nsS1oBNmzZR8nj1ZHHb99bSmL7nqSB278FS888goDApaxtCnI6+Ce45QUFTNj8FByU1S1qzuhQ8nb\nT7cSiIJ+eSx8+FZyi3Jw93Hj8rjoN6yYrz3+ZRzO5MyWjIDBhmVbeGnRq7z70rrmKY3dYefbT9/F\ngJH9cXmcuPu4ySvKYeGvbok5grjghnPiXqO+uoFnHnghSiSq4y6B2ptDwptY8Y93efyeP/PR2l0c\nLjvK8udW8+BnHmFAwBkz0EtjET6i8NhLtEi0QbdMe2+aivLdR3A4HfQfUZyUbEYA9TWN/PK2RVQd\nr8XX4MOd5cLhdnLPH++MKI5zdH8Fhj/AoNEDscUpdx8Mmnzzoh/irY+dadqd5eKLP/0sk+acCtZ6\n4eH/sPwfawgGIqcmLet/+r0Bvn3J/+BrjBxt2J12LrxxDtfevaB5WzILCPdUvMHSHmeb6NVp7202\nYUjJIAaM7Jc0cQBY/LvXqDh0snnU4Gv001DdwF9/9M+I4wYML2bImEFxxQEsu8Xn7rsOuzN2wRpf\no593F6+P2HbJ5+eSnZcVYWdxeZxc/42rIooDl398FIlx7WAgyPY1OyO26ZFE2+iRRHy61SpGqtm4\ndEuUm7RSir1bD1gjij6J+RA0MX3+mZRt2sOKf8Q2Sra8Vl5RLvc9+3XeemYlH60tDS2fns+YqaMi\njsstzI7rzp3fL3q605bnpaZJJHR18pZogQij1aLaHay4fc1/XcHalzdGLc+6s1zMvHxa1PG5fXO4\n+q7LufquUwlr62saWfn8Gj5au4vCAQVc+Nk5jJw4nI8/3BchFC6PK8ofowktEm2jXbaj6ZZTjFRx\n1uVTo5ZRbXYbY6aNwh1yp24vLo+TW//n07jcThwuOyKWOIybVcKUi9qOiag9WccDNz7Mq0++SdnG\nPax7fRMPf+lxZl4xjdFnjsDpdjbXCr3m7gWMmxU/DkRPNxJHL4tadEsjZarwNfh4+MuPc3RfBYbf\nwOF2kJXj4Z4n76RwYOfiDyqPVLHutQ+or25k4rljOX366ITsJy8++gpvP7s6ynDZJy+Lh974PtUV\nNdRW1jNo9IDmmJC20IbLxOmu4eQZFe7dU3D3cfOtp7/KrnVlHNxVTtGQQiafPx67I35l7EQpHFjA\npbde2O7ztryzPUocwFolObL3OEPGDKRwYN8YZ8ZHTzcSp7e7bOspRgtOHq3i5LFq+g0rZtKccUkR\nh86QnR/bxdo0gmTnZXW4XT3dSJze7GClRxAhlFL869evsuIfq7HZLVuB3WHnq7+9HUGorqhm2Lih\nbbp0B/wGuzZ8jGkEOX366HavfLTkos+ex6HS8gifB7vdzvDxQyno37kqWnokkTi91YCpBSLER+/u\nYuXzawj4DAgLy/rFLYuwO+3Y7DaMQJDzr5vNdf99ZUz7wa4NH/P4N55u9pA0gyY33/8pps8/s8P9\nmnrxJA6WlfPmX1bgcDowDZMBI/vxpYdu6nCb4WiRaD+9adqhBSLEyhfXRnkmgjXXDy+Zt/rf7zNs\n/FBmLZgacZy33sfvv/5U1HLmn3/0T0ZMHEbxkMT/kQK+ABuXbaHsgz0UDy1k7qfOZsrciXy8eR+j\nzhjGiAnD2m6kHWiRaD+9RSS0QIRoGSQVD1+jn+V/XxklEJtXxK40aJom7y3ZGDO0OxYNtY384tZF\nVB2rwdfox+Fy8J/fvY7YbThdDpSp+OTdC7jg+rMTai9RwkWiCS0WrdMbwsk7ZaQUkV+IyA4R+VBE\n/iUimZuLrA1mzJ8SUfWqNRpqo+teeOu9zUV0wgkGgjTWNibcj9f++BYnyk82j2YMv4FpKoKBIN56\nH75GP/96ZAk73i9LuM1EiVfRSxOfnh5O3tlVjKXAJKXUZGAX8J3Odyk9zFwwleEThjYbFeOtXtid\ndibPjc6qFJ4RKhx3loszzmu9elY4G5d+iOFvvVq23+vnzWdWJtxmogSNIB+8uYXFv3+DE+/swe83\ntEi0g54oEp2aYiil3gh7uRaIXy0mw3E4Hdz9+y+x5Z2P2LpyOzl9cyga3JcXHn4ZIxDEDJo43U6y\nC/ow/5a5Uef3H1bMBTeczcrnT9kyLI/J0zl9xmlRx8ftR4yEOLGoOhadH6IzWJGsv6XqeE0oktWN\ny+Pgut/fDOjpRqL0NNtEMm0QXwCei7dTRO4A7gAYPCiqtm9GYLfbmHLhJKZceMoF+rQpo1j+3GpO\nHK5k3KwSzv3kzOaktC259u4rmHD2WN5dvI5gwGDGpVM5c+7EdkWcnnvNbF557I2YBXmacLjsSU+Q\n89Ki16g4XNnslOVr9OH3+Vn189e56pc3aANmO+hJtok2Xa1FZBkQK9njfUqpl0LH3AfMAK5VCfhu\nZ6qrdSYQNII88a2/sDNkYzBNhREwaMo+63Da6ZPXh/ue/Rq5fZOXXfqbF/2I+pqGqO12u51frvgR\n+/1WMhstEu0nHe7aXeZqrZRq1fwuIrcCVwIXJyIOmtaxO+x85Ve3cmDnYfZ/dJCCULaqt/+2iuqK\nGiaeO46LbzovqeIAWLn74u1CL4V2hu487ejUFENELgO+BVyglIr++dF0mGFjBzNs7KlCOhPPSW3O\nzbMun8KqF9dGGEhtdhtjpo9qDgLTItFxuuu0o7OrGIuAXGCpiGwSkceS0Kcex+HdR3n9qbd585mV\nMZPPZgJXLZzPgJH9cfdxY7fbcWe7ySvK5eb7PxVxnF4C7TjdMR+mDvdOMf/6zauseHY1QcNE7IKI\n8Nn7rmXWguhkMenGNBXb1+7iUOkRiocUMvmC8XGTAeuQ8eSQqnyYvTonZXdhz9YDrHjuXfy+AMFg\nEMNvEPAF+NuDL1BflXkzMptNmHjOWObfcgHT5p3RaqZwPZJIDpk+ktACkUI2vLGJgD96udJmt7N1\n1Y409Ci5aJFIDk0ikYnh5L0yFqO+uoF3X1rHvo8OMnjMAOZcM4u8ohRV5urh6zracJkcMjWcvNeN\nIE4cPsmPrvslLz++lI3LPuT1P73ND6/9BYfKjiT9WjPmT8Hpjk4DZwaDTJozLunXSxfhI4mmh6bj\nZJLLdq8TiOcffpmGmobmCuEBn4G33sfff/Ji0q81ctIwLrzxHFxuJ3aHHafbgdPt4LPfu77VYrzd\nkfBAr5K8Qi0SnSRTRKLXTTG2r9mJaUaP+/du3Y8RMJJWwq+Jq++6nJkLpvHhO9txuh1Mu/iMqExQ\n1RW1vP7U22xfs5Ocwhzm3XQBZ14QHRDWnSjJK6RUTzs6RSb4TvQ6gXC6nfi90YZDsduQDta+aItB\nowcwaPSAmPtqKuv4yWcfoaGmkaAR5Oi+Cg7uOMxlX7yoQ0luMwktEp0n3baJXjfFOPsTZ+F0R+qi\nw2ln+rzJ2O1d/3G89beVNNY2RhTA8TX6efUPb8at69md0NON5JGOaUevE4irFl7CmKmW+7A7VHBm\nyOmD+fS3r05Lf3asLcWIkdbe5rClxHCaDrRIJI+uFoleOcX46qLbOVRazuHdR+g/rJjhE4YmtQhw\ne+g7sID9Ow5FbQ8aJvmpWnpNA3q6kTy6Mvir140gmhhSMoizLpvKiInD0iYOAPNuOi8q1Z3daWfE\nhKEUD03/Ongy0SOJ5NFVtTp6rUBkCqdNGcWn772arBwP7mw3DreDMVNG8eVffD7dXUsJWiSSR1fk\nw9TBWhmCETA4uu842fnZFPTLS3d3Uk5pTaWebiSZ8MQ0Olirh+FwOhgyZlCvEAfQI4lUkIqRhBYI\nTdrQIpF8wv0mkoEWCE1a0SKRfJIpElogNGlHi0TmogVCkxFokchMtEBoMgYtEplHr/Ok1GQ2TR6X\n4ejl0PSRlBGEiHxDRJSIFCejPU3vRhcRzhw6LRAiMgyYD+zvfHc0mki0SKSXZIwgHsYqntPDsy9q\n0oUWifTRKYEQkauBQ0qpzQkce4eIrBeR9ZWVJzpzWU0vRItEemhTIERkmYhsjfG4GvgucH8iF1JK\nPaGUmqGUmlFYWNTZfmt6IVokup4OF+8VkTOAUcDmULj0UGCjiMxUSvWMTCeajEOn2e9aOjzFUEpt\nUUr1V0qNVEqNBA4C07Q4aFKNHkl0HdpRStMt0SLRNSRNIEIjiYpktafRtIUWidSjRxCabo0WidSi\nBULT7dEikTq0QGh6BFokUoMO1tL0GMKXQJvQS6GdQ48gND0KHeiVXLRAaHosWiQ6jxYITY9Gi0Tn\n0AKh6fFokeg4WiA0vQItEh0jLZW1RKQW2NnlF45PMZBJXqCZ1h/IvD7p/rTOWKVUp6s/p2uZc2cy\nyoIlCxFZr/vTOpnWJ92f1hGR9cloR08xNBpNXLRAaDSauKRLIJ5I03XjofvTNpnWJ92f1klKf9Ji\npNRoNN0DPcXQaDRx0QKh0Wji0iUCISI/FJFDIrIp9FgQ57jLRGSniJSJyL0p7M8vRGSHiHwoIv8S\nkYI4x+0VkS2hPidl2ahF+62+X7H4dWj/hyIyLdl9CLvWMBF5W0Q+EpFtInJ3jGPmikh12PeYUEbz\nTvar1e+giz+jsWHvfZOI1IjI11ock9LPSET+KCLHRGRr2LZCEVkqIqWhv33jnNv++0splfIH8EPg\nnjaOsQO7gdGAC9gMTEhRf+YDjtDzh4CH4hy3FyhOUR/afL/AAuBVQIDZwHsp/I4GYSUdBsgFdsXo\nz1zg5a74n0n0O+jKzyjG93cEGNGVnxFwPjAN2Bq27efAvaHn98b6f+7o/ZVJU4yZQJlS6mOllB94\nFrg6FRdSSr2hlDJCL9dipezvahJ5v1cDf1YWa4ECERmUis4opcqVUhtDz2uB7cCQVFwryXTZZ9SC\ni4HdSql9XXCtZpRS7wCVLTZfDTwdev408MkYp3bo/upKgfhqaAj4xzhDoCHAgbDXB+maf9AvYP0C\nxUIBy0Rkg4jckeTrJvJ+0/KZiMhIYCrwXozd54S+x1dFZGKq+0Lb30G6/m9uBP4eZ19Xf0YDlFLl\noedHgAExjunQ55Q0V2sRWQYMjLHrPuD3wI+xvuwfA/+LdWOmjNb6o5R6KXTMfYABPBOnmTlKqUMi\n0h9YKiI7QgreYxGRHOAF4GtKqZoWuzcCw5VSdSE70r+BkhR3KeO+AxFxAZ8AvhNjdzo+o2aUUkpE\nkua7kDSBUHEqcLVERP4AvBxj1yFgWNjroaFtKemPiNwKXAlcrEKTtBhtHAr9PSYi/8IapiXrnzOR\n95vUz6QtRMSJJQ7PKKVebLk/XDCUUktE5HciUqxSWO4gge+gSz+jEJcDG5VSR1vuSMdnBBwVkUFK\nqfLQ9CpWyGqHPqeuWsUInxNeA2yNcdg6oERERoUU+kZgcYr6cxlWRfJPKKUa4hyTLSK5Tc+xDJux\n+t1REnm/i4HPhyz1s4HqsKFkUhERAZ4EtiulfhXnmIGh4xCRmVj/PymrxJzgd9Bln1EYnyHO9KKr\nP6MQi4FbQs9vAV6KcUzH7q8usvj+BdgCfBjq1KDQ9sHAkrDjFmBZz3djTQVS1Z8yrPnYptDjsZb9\nwbL2bg49tqWiP7HeL7AQWBh6LsBvQ/u3ADNS+JnMwZoCfhj2uSxo0Z+7Qp/FZizj7jkp/r+J+R2k\n6zMKXS8b64bPD9vWZZ8RljCVAwEsO8IXgSLgTaAUWAYUtvx/jvf/1tZDu1prNJq4ZNIyp0ajyTC0\nQGg0mrhogdBoNHHRAqHRaOKiBUKj0cRFC4RGo4mLFgiNRhOX/w/UZO0BYb+zbAAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f67c6ef8128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_train, y_train = create_toy_data(add_class=True)\n",
    "x1, x2 = np.meshgrid(np.linspace(-5, 10, 100), np.linspace(-5, 10, 100))\n",
    "x = np.array([x1, x2]).reshape(2, -1).T\n",
    "\n",
    "feature = PolynomialFeatures(1)\n",
    "X_train = feature.transform(x_train)\n",
    "X = feature.transform(x)\n",
    "\n",
    "model = SoftmaxRegressor()\n",
    "model.fit(X_train, y_train, max_iter=1000, learning_rate=0.01)\n",
    "y = model.classify(X)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)\n",
    "plt.contourf(x1, x2, y.reshape(100, 100), alpha=0.2, levels=np.array([0., 0.5, 1.5, 2.]))\n",
    "plt.xlim(-5, 10)\n",
    "plt.ylim(-5, 10)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3.5 prohit回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一般线性模型框架:$$p(t=1|a)=f(a)$$\n",
    "其中$a=w^T\\phi$\n",
    "对于每个输入$\\phi_n$,计算$a_n$,然后设定目标值$$\\begin{cases}t_n=1,&\\mbox{if}  a_n\\ge\\theta\\\\t_n=0,&\\mbox{other}\\end{cases}$$\n",
    "如果$\\theta$从概率密度中抽取，那么对应的激活函数由累积分布函数给出$$f(a)=\\int_{-\\infty}^ap(\\theta)d\\theta$$\n",
    "\n",
    "作为⼀个具体的例⼦，假设概率密度p(θ)是零均值、单位⽅差的⾼斯概率密度。对应的累积分布函数为\n",
    "$$\\begin{align}\\Phi(a)&=\\int_{-\\infty}^a\\mathcal N(\\theta|0,1)d\\theta\\\\&=\\frac{1}{2}\\bigg\\{1+erf(\\frac{a}{\\sqrt2})\\bigg\\}\\\\\\end{align}$$\n",
    "\n",
    "其中误差函数$$erf(a)=\\frac{2}{\\sqrt\\pi}\\int_0^aexp(-\\theta^2)d\\theta$$\n",
    "在实际\n",
    "应⽤中，使⽤probit回归得到的结果倾向于与logistic回归得到的结果类似"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 Laplace 近似"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "后验概率不再是高斯分布时，需要以某种方式来近似,先考虑单一变量z，分布$p(z)$定义为：$p(z)=\\frac{1}{Z}f(z)$,其中$Z=\\int f(z)dz$\n",
    "\n",
    "目标是找到一个高斯分布的近似$q(z)$,其中心为$p(z)$的众数。首先，寻找$f(z)$的驻点$$\\frac{df(z)}{dz}\\bigg |_{z=z_0}=0$$\n",
    "然后泰勒展开$$lnf(z)\\simeq lnf(z_0)-\\frac{1}{2}A(z-z_0)^2$$\n",
    "其中$$A=-\\frac{d^2}{dz^2}lnf(z_0)\\bigg |_{z=z_0}$$\n",
    "在$z_0$处有概率分布的局部最大值，所以$f(z)\\simeq f(z_0)exp{\\frac{A}{2}(z-z_0)^2}$\n",
    "指数部分为z的平方函数，所以可以使用归一化的高斯分布近似得到\n",
    "$$q(z)\\simeq (\\frac{A}{2\\pi})^\\frac{1}{2}exp{\\frac{A}{2}(z-z_0)^2}$$\n",
    "对于高维空间有$$q(z)=\\frac{|A|^{\\frac{1}{2}}}{(2\\pi)^{\\frac{M}{2}}}exp{-\\frac{1}{2}(z-z_0)^TA(z-z_0)}=N(z|z_0,A^{-1})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXZ/YkLGENq4CAIoioRUTEFXdQqm39SdW6\ntLXeaq9dXKq13l57e2977XattZZaa+tSW3fcl7orIKCIbAoBVJB9J8usn98fZxImySSZZCaZOZPP\n8/HgIXPmzDnfQfLm+z3fTVQVY4wpJp58F8AYY3LNgs0YU3Qs2IwxRceCzRhTdCzYjDFFx4LNGFN0\nLNiMMR1ORO4RkS0isrSZ90VEbheR1SKyRESOTHnvDBH5KPneDzO5nwWbMaYz3Auc0cL7ZwKjk7+u\nAP4AICJe4PfJ98cCs0RkbGs3s2AzxnQ4VX0D2NHCKTOBv6ljHlAuIgOBScBqVV2jqhHgoeS5LfLl\notBt1advLz1g2OB83NqYLmHxe8u2qWq/bK4x7bSpumP7rkzvtwyoTTk0W1Vnt+F2g4HPUl6vTx5L\nd/zo1i6Wl2A7YNhgXnnnn/m4tTFdQu/QuE+yvcaO7bsy/jntHRpXq6oTs71nruQl2IwxppENwNCU\n10OSx/zNHG+RPWMzxhSCOcDXkr2jk4HdqroRWACMFpERIhIALkie2yKrsRljOpyI/B04EegrIuuB\n/8CpjaGqdwHPAmcBq4Fq4LLkezERuRp4AfAC96jqstbuZ8FmjOlwqjqrlfcVuKqZ957FCb6MWVPU\nGFN0LNiMMUXHgs0YU3Qs2IwxRceCzRhTdCzYjDFFx4LNGFN0LNiMMUXHgs0YU3TyEmzheJQ1e7bm\n49bGmC4gL8EmEgKwcDPGdIi8NUVDvtGAhZsxJvfy+ozNws0Y0xFytrpHctOFhcAGVZ2R6ef2h9sq\nAA7skdVqxsYYk9Ma2zXAivZ+2GpvxphcyUmwicgQYDpwdzbXsXAzxuRCrmpsvwWuBxLNnSAiV4jI\nQhFZuHPb9mYvZOFmjMlW1sEmIjOALaq6qKXzVHW2qk5U1Ym9+vZp8Zoh32hCvtGs2bPVAs4Y02a5\nqLEdC5wjIutwNjM9WUTuz8F1rfZmjGmXrINNVW9U1SGqOhxnB5lXVPWirEuWZOFmjGkrV8wVtXAz\nxrRFTnepUtXXgNdyec06Nt7NGJMpV9TYUlntzRjTGtcFG1i4GWNa5spgAws3Y0zzXL0TvD13M8ak\n49oaWyqrvRljUhVFsIGFmzFmv6IJNrBwM8Y4XP2MLR177maMKaoaWyqrvRnTdRVtsIGFmzFdVVEH\nGzQMNws4Y7qGog822L++G1jtzZiuoEsEWx0LN2O6hi4VbGDhZkxXUHTDPTJhQ0KMKW5drsaWympv\nxUE1gYbnoft+h+67F42tz3eRTJ51yRpbqpBvNLWxVazZs9Vqbi6kGoPdt0J8NWgt4IHa59Fu/4aE\nTsh38UyedOkaWx3bFcvFwm9CbFUy1MDZATIC++5CE7UtfdIUMQu2FNY0daHwW0C46XHxQGxFpxfH\nFAYLtkYs3FxGQi28F+y8cpiCYsGWhs1WcJHQaUCaAJMA+MZ0enFMYbBga4bNVsgtja1H9/wc3X4J\nuvM7aO2rqGrW15XABCiZAfic2puUgHSDHjcjYn+9u6ou3yvaGus1zZ7GNsHuG0BrnAPxvbBvNsQ3\nQdmsrK8vZReiodMh+qETaoHDEfFnfV3jXvZPWgas5pal2kdAGz/gD0PNk2iiJie3EG9fJHQSEjzK\nQs1YsGXKhoRkIfoRzjCMRsQLiU2dXhyTHyJyhoh8JCKrReSHad6/TkQWJ38tFZG4iPROvrdORD5M\nvrewtXtZsLWR1d7awTMw/XGNgvP31hQ5EfECvwfOBMYCs0RkbOo5qnqbqh6uqocDNwKvq+qOlFNO\nSr4/sbX75eUZWzgao3LzDkZWuPMvtT13a6PSL8HuJUAk5aAfAhMRb8+c3UYTeyGxC7wDrDmaA+F4\nNJf/gE8CVqvqGgAReQiYCSxv5vxZwN/be7O8BFvQ49y2crMTxm4MOJtInznxH4x2+x5Uz4bEXkAg\ndDyUfTMn11cNw97fQeRdSAaals5CSmbk5PpdlUio/u95Bvo2aiLOVtXZKa8HA5+lvF4PHJ3+vlIK\nnAFcnXJYgZdFJA78sdG1m8hbr+jIMicIKqu2Wu2tC5DQ0WhwEiR2g6cUkUDuLr73/5xQIwEac45V\nPYB6+iPBSbm7j2nJtkyaiBk6G3i7UTN0qqpuEJH+wEsislJV32juAnl/xjayrB8jy/pRuXlHfQ3O\nbey5W2ZEBPGW5zTUNLYWIvNo2jkRhppHc3Yfk7UNwNCU10OSx9K5gEbNUFXdkPzvFuBxnKZts/Ie\nbHXqa3AuDTjrNe1YGluH7v0VuvMadO9v0ViyVbPv3uY/lNjeKWUzGVkAjBaREeL8y3YBMKfxSSLS\nEzgBeDLlWJmIdK/7PXAasLSlmxXUAN3GzVNw3/M3a5rmnkaXw+6f4nQ+KMTXQ2Q+2uPWlie6+8Y2\n/57pVKoaE5GrgRcAL3CPqi4TkSuT79+VPPVc4EVVrUr5eAXwuIiAk1kPqurzLd2voIKtjtsDLjXc\nwDoWslZ1Nw1X8FBnwG/VPc6c0Lrnao2VZj+rweSOqj4LPNvo2F2NXt8L3Nvo2BpgQlvuVTBN0XTq\nnr8Brmue2lzT3FBViK1L/2ZsNQRPBhoP7RAITEZ8zYyfM0Uv62ATkaEi8qqILBeRZSJyTS4KlsrN\nHQwWbtkREWdie9o3y6DsIvCPAQLJCfBB8I2Gbt/p1HKawpKLpmgM+IGqvpd8wLdIRF5S1eYG3rXb\n/tqbExJuaZ7mqmmqsU1QfV9ysneps6pFaDrJZw/FK3gm1D5NwwG+QSiZ4fSw9vxPp3c0th58gxDf\nyHyV1BSIrGtsqrpRVd9L/n4vsAJnMF6HcWMParZNU43vhN3XOUMbdB8ktkDV/VB1V+sfdruyWRA6\nDvAna28+CJ0EJefVnyK+EUjoOAs1A+S480BEhgNHAPNzed103NrB0O5e09qnkytkpK5hFoHaV9GS\nCxBvr1wXtWCIeKHbVWjpxU6geyoQT/d8F8sUsJwFm4h0Ax4Fvquqe9K8fwVwBcDAwUNydVtXBly7\npmNFV+K0+huRAMQ/hSIOtjri6QGeHvkuhnGBnPSKijPj+FHgAVV9LN05qjpbVSeq6sRevXMfPG7s\nQW1T09Q3mLT/uzQGnv65LVgn0to3nBV1t1+I7roZja7Kd5FMEchFr6gAfwZWqOqvsy9SdtzWg5rx\n/gqhs2lawfaB/2DXDmvQmqdh3x8gvsFZXTe2HHb/GI1V5rtoxuVyUWM7FrgYODllkbizcnDdrLip\ngyGTjgXxDYUePwLvAJyA80HwGOjeZL0+V1CNQfXfabp1XgSqHsxHkUwRyfoZm6q+BRTkeAO3PX9r\nrWNBAuNR/++dXlEJuXvNscQu0Hj69+JrO7cspugU5JSqXHNTwLXWseAMWC2MHsF4LM6axZ8QTyQY\nOWEY/mAbgralTgBP505BU01AZAFE3nGGkwSnIf6M1yEzBahLBFsdN60BV+iT6Ve9t5bZ1/2NeKxu\nuSDlsv+axfjjDsno8yIBZ2ep2hdoOPA2AKUX5Lq4aWlsI0QWQeRfENuE0ywWqH3NWaiydGanlMPk\nXkHPFe0obulgKNSlkGr21fKH7/6Fqt3V1FbVJn+F+fOND7Bry+7ML1T2NQhNT+7Y7gNPL+h2FRI4\nosPKXkerHoBd34Xqv0LsE/Y/61MgAtUPovFdHV4O0zG6ZLDVcVPAQeHMN1386lKUppsdJxLKghcW\nZ3wdES/S7WLofT/0vgd63Y2EjstlUdPS6AqoeRqIAs0858PrTF0zrtSlg62OG3pQMx4W0glq99am\nNEH3i0ViVO2qbvP1RLyIp1vnzXkNv0bT3thGRJI1SeNGFmxJjQf4FmLAFcpSSAdPGoUnTQgFSwKM\nm3JwHkrUVs3V0hrphCax6RgWbI24JeAgf7W3QaMGMGn6kQRL9u9dECwJMGbSaEYdOaLTy9NmgalA\nc7WxEmfllB43u3s4TRfXpXpF26LQe1Drwi1fPaezbjyXQ48dwztzFpCIJTh6+pEcMW28O5ZQ8k+A\n0BSofQenSZr8MQidCv4jIDDBQs3lLNhaUehrwOVrGXIR4bATxnLYCe7bV0BE0LKrIXiqM35NQhA4\nHvFV5LtoJkcs2DI0sqxfwQ7wzXftzY1ExFl51z8m30UxHcCesbVBoT9/y/ezN2MKhQVbOxRywBVK\nz6kx+WRN0SwUYgeDohB+k0Dt86C1rA4fjASOZWTP9i/uqRqG6HLAA/6x9mDdFDwLthwopA4G3fdH\nCL9O3QDUUHwDtZH3qZRrGNmj7eu2ae18qPo/9i/gImj3G8A7GBLbwDsY8ZTlrPzG5IIFWw7lu4NB\n45sg/CrOVKE6UUK6B4lvZs0e5393pp0LGt8O+35Dw0nqwJ6fAF4QP2gMDc2AsovcMdTDdAl5CbZo\nJMZna7cAMHSEe5e1TievSyTFVoJ4QaMNj2stWvsiAe8QaklQmTga8ZS3HnDhN4CmU6ecieKx/Tuw\n1z4L3gooOS0HX8KY7OWl8yDo8zGyr/PD/tnaLfW/ikleOhg8LWzoElsF4ZcJhd8guPdXaGxl650L\nWkXaDWSaCEPtk20pqTEdKq+9oiP79q7/BVjAZcs/3tkdPa26+ZExIEKw+mGURMtDQwJHZj4RPLEX\nTdSg4blo7VtoYl8bC29M7hTMM7a6cAOoTAm3YmmqdkYPquCBHreie/8H4ltAPKAR0jcnY4TwIC3N\nXPAdAv4jIfIera6G4R0IOy537gmwL4Z2+zYSOiEXX61ZqmFnj4TwK04TPPAFKLsU8dog5a6sYIIt\nVV3IVW7b0aAGVwwh19E9qOIdiJTfjsY/Bw2jVX92dn9qQgFnEntzMxdEBO1+rbP7fPgNwAO+YVD9\nGPs7FLzOdWLrgGjD/Zz33Yn6x3ZsyOz+GcQ+or7DJDIPosvQXncintKOu68paAUZbHUa1OKKLOQ6\nugdVvIOc34ROR/dV0qTG5SkH3wENDqWbdyoizm5YwWPqz9PAF5xwS3wOvjHg6QvVj6QphTOmjtLz\ncvjNUq4eq3SeHTboBVbQsFODK5nRIfc1ha+ggy1VcyHn5oDrlB7U4FRnJdjw64DHWUCRANLtJiTN\n5mKNa28AI8qA2pcgsRX8h0FwKtLjuvrPaM1TpF/jLAZam5OvofEdUPs0RFc4m0eHZkLsU+f7NFnM\nNwyx1Tm5r3En1wRbqmJrqnZkwAmCdPs2WvJFZ/aApycEDkdoefZAXcDV1D5H5ed/YWTJLiAGkflQ\n8xja838RT4lzsv8LwANpruIB30FZfweNbYLd1zk1MWJOLa32bSi7qJlPBMA7LOv7Gvdy9VzRYutV\nbdyDmkviHYSETkECR7UaanUUJVjzMCFiVNaUUlnTw6mBxTdDzf7hHeIb5GzK0kQCqv6AJmqyK3z1\nfaDV7B96kgDCzr4F3oE0+fdZfBCalt09jau5OthSpYac28fGFcwmM/ENoE4ohZJ/VSprelBZUwqR\ntxueG5xE2gZAotp5zpaN6BLStDchsR26XZd8/ucDxOnJ7fk/SEv7lpqi58qmaGuKpama7ylaSBB0\n/1CRunCrJUFldQ883pS13+KfkP7fyTDEK7MsR1lysHAa3j5I9++h3a4B4jZB3wBFVGNLpxiaqvlc\nIkk8/cA7BBp1MoQooaTkXCBlaSTPwP1j2BoIJK+RhZIZ1A1N2c8HwWPqg0zEY6Fm6hVljS2dQqvF\nbduwgw2rNtJ3SG8Gj2p91Y18zUGV7teje34MmpxJoHGn6ReaRigZeGv2rEIZwIGefhD/nIY9pD7A\ni+661nn4H5gKJefs73jIRGg6xNdD7SsgAWcgrn8slF2Zq69pikyXCbY6+R42Eo/FufeWh/jw9eV4\n/V7icWXoQQO56v8uI9St9R/2zl4DTrwV0OsuZ8hIYgf4xiDehkFcP/7N8z00/CgjffOARHKMW3eo\nuo/6cXQ1j0HkbbT8VxnXsEQEul2Jlvw/J+A8/RDfgNx+UVNUulywpcpHLe6lv73Oh2+sIBqJEY04\nvXyfrljPg//zOJf/7KsZX6cz14ATPM7OTi2oH/8mZazhIiDBgaXAru/SeBklEtsg/BaETmpbOby9\nwNvCRH9jkrp0sNXpzFrcG4/MJRpuuKxQLBrng1eXEovG8Pnb9r+kkBa5hEYDfHctRGvLGRlqNMle\nayH6QZuDzZhMFXXnQXuk62zIZYdDpCaS9ngiocSiGe5Qnkah7cEQ8o0m5He25nOGiKQOv/CBxyap\nm45jwdaMjupRHTP5IMTTdCpTxfB+hEozXCKoGQW3yYx/PCHpSQivU6b6gPM6e3oa00Es2DKQy4A7\n75qzKO1egj/gNDm9Pi+B0gAX3vzlnJQV2hZwmz/dyvN/eYVn//wvPl+zOWdlAOfZnPT4KXiHESJI\nSEpBurPG8x3E564xhcZd7BlbG+RizbjeA3pxyyPX8uaj81i75BMGHFjBCecfQ5+BuX821loP6ssP\nvMlTf3iBeCwBKC/+5VVOvfREpn/jlJyVQbwVSPmv0fgWIEyJdwiCsGbPKqDzdq43XUtOgk1EzgD+\nD2dxrrtV9ee5uG4hS9ejmmnAdSsv48yvd95cxnQdDNs27uSpO5+v75kFiMQTvHTvaxw5bTwDR1Tk\ntAzibfhns3/vUwu4rqK1nBCRE4EngbXJQ4+p6q2ZfLaxrINNRLzA74FTgfXAAhGZo6rpVjcsOoU2\n8LclqVO0Fr68OO05sWicxa8uyyjYNPYJWvWn5CYyIQiehpTNyniSPVjAdRVtyIk3VXVGOz9bLxc1\ntknAalVdkyzEQ8BMoEsEW518D/zNVF3tbRFCtDxELBbHtyNl9Q2P4EnTudGYJrahe27cv96aVkPt\ns2hiE9L9+jaXywKu6GWTE23+bC6CbTDwWcrr9cDRjU8SkSuAKwAGDRycg9sWrmyaqZ3ltGOP4p1f\nvkpsgIdYb2fGg29HDV6PcPi08a1+Xmue3r/9Xr0IRBah8c3OjIV2aBxwChwY2uAsHOnpA8EpSKYb\nzJishKOxtvSs9xWRhSmvZ6vq7JTXGeUEMEVElgAbgGtVdVkbPluv0zoPkl9yNsD4QyekWYOm+BRy\nwJX378mXrzuHR26bAyJEBvugfzdmXnQiFUP7tn6BWCVpt+YTPyQ2OPuMtkITu9GaJyC6CKQnUnIO\nEjgKcAJOCVOz6zoq92wGjTCyNAxV96I9f4b4spxYb1oV9Pjqa/gZ2KaqE7O85XvAAaq6T0TOAp4A\nRrfnQrkItg3A0JTXQ5LHTFKh7sA19YtHc+iUMSx+bRkklMNOHMvO7rHM5qD6DnSerTVeElyj4BnU\n6r1V96C7vweJfTgBuR7dtxotOR9PcuUQrX6CUHwzEKUWpbLaWeFjpPc3UP6rNn9fk1et5oSq7kn5\n/bMicqeI9M3ks43lItgWAKNFZETyZhcAmU967GIKrRZX3r8nJ54/pf51XZy1NkVLSmag4Zec1T6S\nEgkf+6pGUta9D35vy/fVmmdSQq3uYBiqH0JDZyBSAuFXqZtn2mAtuL27ENYysnxEm76ryatWc0JE\nBgCbVVVFZBLOONvtwK7WPttY1gN0VTUGXA28AKwA/plsF5sWFPo6ca0N8BVPP6THz8B3MKoQrhHe\neqYH//1NLzec9lNWvb+2yWcaiLxP+qasD+J1n236xCKEJzmTQVve7NkUlOZyQkSuFJG69ae+DCwV\nkQ+A24EL1NHmjBHVzn/cNf7QCfrYP5/t9PsWsspt+8Mjmxrclk+3UvnBJ3Tv3Y2xkw/C481+ckll\n1f7waFyD27B6I7ddegeRcMOQCpYG+fkLPyYYSj/0I7H3f509QJsIIOW/QbwDSVT/HWqeoOHqIIB3\nOJ7yXwPOZPs61pO6X+/QuEXZPvMad9gE/eczL2Z07qEHDMj6frlkMw8KRLZNVNUED/zXoyx84QPE\nI3hECJQE+O4fr6BiWHZN3ZYWuZz71EJi0aY7zQuw/J2VHHFy+h5WKTkbjbxPw/1OveAbXr/em5Se\n55yTWJ8cVhIC8SHdv1f/ibpeVNjfkwoWcl2dBVuBaW/AvfvcYha9tKTBkkjhmgh/vPY+bnn4+zRe\n3rtt4mjNoxwYfhq0msrYOCrjMxk5aAK11RESiabBpkC4Jtr0UkniOwQtuwKq73bO1jj4RiM99o+B\nE4JQ/guILnaGe0jf5HLgobTXTLcnqgVc12TBVqDS9aS2FHBvPjq3yZJIqsqOTbvY/MnWrGptuu+u\n5E5TTu1qpG8JxD6mcuPNVEw6kMCLHxCpbrjTfDwWZ8zRLffUe0InoaGpEFsPnu6IxxlmsmvbXnZs\n3EH/Yf3o1qMU/Ec4vzJkAWcs2FygvhbXQsBFa9M8iAc8HoiG07+XEd2T3EW+Ue1Loxwor6Djz2fw\nlFGsW/YZunEP4hH8AT9nfXMa5X27t3p5wQ8+p3czGoly7y3/4MM3V+AP+IhFYkw9bzJf+v4MPNL2\nGqcFXNdlweYiLQXcxNMnsPmTrU1W5/UH/QwalcX+APGN+zdQafgGGlvDqPIKrr3pEl58ewEr311N\nIODjzLOPZvi4oWkv15KHf/UUS99aSSwSI5acnP/2E+/Sd3AvTrpgaru/gj2H63psPTYXarwxNMAJ\n50+hYlg/AqXOIFav34s/5OfSWy/A48nif7OnArTpqr/xOCz6127mPb0QEQ+nTz2aa75/IadddjLx\nvmVtXuQyFosz/+lFTYI5UhvhXw9kueFyipBvdMq0LRsuUqysxuZSqsqgUCnB0iBrk+F23V+vYvEr\ny1g57yPKK8qZMvMoeg/IcvMTTzkEj4XwXFJ7MGNR4fkHe7J9y5P07NuDQyYfBLR/F61YJEY83rQT\nAqBqT03a49mwZmpxs2BzoTcfm8+cO1+gdl8twZIAZ3z9ZKZdeBxrPttBxeiBTDyt5R2l2kq6XUUk\n2oN41dMEQgnWrw7x8J392PRpAIjwwl9erQ+2Om3dZCZYGqDPwF5sXb+94b0FRh7ecTMMrJlanCzY\nXGb+M+/x6K+fIlLrNNmq99bw9F0v4fF6OHmW8xwqk17UtvGxY/d0/veSNURqwqg2fJC/Y9OuZj+Z\nacAJwqwbz+WuH/yVaCSGJhSP14M/6Oe8a87K0fdomdXiioc9Y3OZp2e/VB9qdSK1EZ778yvUzSLp\niKlafQc712wcauLxcOCEYa1+PpM9GMZMGs2191zFkadMYPBBA5ky8yh+9PfvMujA3K7m25q653Ah\n3+j653D2LM5drMbmMrs27057vHp3FYl4Aq/PmX2e68n2Pr+Pc759Bk/e8TyRWqczQTxCIOTnrG9m\ntkdCSzMY6gwZPZCv/2xWu8uZa1aLcycLNpepGN6Pz1dvanK8V//y+lBLlcuAO/H/HUuvinJe+Msr\n7Nq6h5GHj+DsK0+lfybrt6WWKYOAKzTpAg4s5AqVBZvLnPvvZ/Gn6+4jkjIswh/0MWTMIH488xcE\ngn6O//Jkpn5pMt6UCfC5CrgJJ45jwonjsvwWyTK1swc1n1I7G6wWV7gs2DqBqkJ0GRqZC/iR0AmI\nr309feOmHMwVv7qEJ+94ji2fbKXXwHJq9tSyfO5HxCLO2miP3/4cqxev4+v/3XTJqvytB6eQ2A2e\nEqDh0t5t7UEtFFaLK1wWbO302crP+fCtFQRCAb5w6mH0quiZ9jxVRavugPA7zkKKCFr7PFo6C0/p\nzHbde+zk0Yyd7PxQvTNnIQ/f9mR9qIHTmfDh68vZuHYLA5sJrcYB16HhFl3kzDdNJJ8PBicjZf8G\nUtKwTCm7aIH7Ag6sFlcoLNjaSFX5521zmDtnAbFoHI9XeOoPL3DRLV/mqNMPb/qB2IqUUANn3YsI\nVD+ABqci3j5p75GIK15f653WqxZVEq5pOjNAPMK6pZ82G2x1MpmHmpXYGnTPbTRYnig8D01UIz1u\nbloeFz5/S2Xj4gqDBVsbrVq0lnlPLawfcpGIA8S5/9ZHGDflYEq7N6yFaHhuSqil8kL0PfCeWn8k\nGonxxO+e4+3H3yUajjD4oEFccMMXOfCw5odT9BnUG1/A26DGBoBH6NU/fS0ynZF9e3dI81RrH6fJ\nBHqiEP0QElvBk/4H3u0BB9ajmk82jq2NFr64uH64QyqPz8OKuR83/YAESP/HLECgwZF7b3mItx6b\nT6Q2giqs/+hzfnfV3Wxe1/wYqmO/OAmvt2FvqMfroVuPUkZPHJnBN9qvQ5Yrj20E0kyVEj8ktrVe\nprJ+DcbAuZGNi+t8FmxtJCKkW7RRnDebHg+eSPqKcQKSW80B7Ny8m6VvrmgyCTwaifHy/W80W55e\nFT256vbL6T2gnEDIjy/g44BDhvC9P32rQa9oWzQOuKz4x5D2+2sEvJlvoVcXcC0N8HWDdJPwLeRy\nz5qibTTpzMOZ/8x7TWptiXiCsccc1OR88Q1Fyy6Gqr+BJGtWqkiPaxFPaf15W9dvx+f3NVk7LRFP\nsH7VxhbLNOqIEfz0qR+y/fOd+IN+emawDlomcvH8TUrORcOvgSbYX3MLQeh0kLaX060dDI1Zh0PH\nsmBro5GHj+D4r0zm9X/OdUb6ez0ocMmtF1DSLf2S1Z6S6WhgirPENT4ITEQ8DZ/FVQzrRyzadEFI\nr9fL8LGt12xEpH7aU65l9fzN0wfpeRta/SBEl4CnO4RmIqHMZiukLU8RPH9LZR0OuWfB1g7nXTOd\nY86eyNK3VhIoCXDEyYfSo0/LtQ/x9gLvSc2+37Nvd448dQLvv7ykwVxQX9DHtIuOz1nZ2yur4SHe\ngUj3H+S+TEUWcGBj43LFgq2dBh5YwcAcT86+6OYv03dwb17/xzvUVNUycsJwvvKDs+trYrFojKfu\neok3H5lLuCbC8HEHcP715zDskMyfVTVHVVn13lpWzP2Y0p4lHHX64ZSn6VXt8OEh7VDMAQfWVG0P\n21fURf5804MseX15gw6GYGmQmx68hn5Dmo6Hy1Q8nuDuHz7AynkfE66J4At48YiHr//8QsYfd0iz\nn8vVXqi0SA3yAAASB0lEQVS5VrcPqtvDLZ1M91Ht6vuKWq+oS+zcvJslry9r0msaC8d4+b7me00z\nsfiVpfWhBhCLxImEo9xz89+b3C9VTntPc6hYelDTsV7VzFiwucSWz7bh8zd9chCPx/nso8+zuvb8\nZxaln70ArH5/Xauf74j133IhkzXg3Cp1bBxYyDVmweYSFQf0bbbX9IAxg7K6drrljup4MpjWBYVf\ne4PiDDhIH3JdnXUeuER5/54cduKhLHmtYXO0tV7Tzys38/jvnqXy/XWUlZdyysXHc/yXJicHGjsm\nnz2RlfNXNam1ebweRh0+vE3lLMTOBSjODoZ0UjsdujKrsbnIpf95PifNOpaSbiE8Xg8HHjaM783+\nVrMdB1s+28Ztl/+e5W+vpLaqlu0bdvD4b5/liTueb3DeYccfwsQzjyAQ9OMLeAmWBAmWBvnWLy9p\nsTbXkkKsvUFxTNEyrbMam4t4fV6+ePWZfPHqMzM6/8V7XyVaGyW14ztSG+G1h97ijMtOpKSbM0hY\nRLjwpvM48fwpfPTuakq6hzj8pHH177dXodbewL1rwJnMWLC5XDwWJxaJESwNNnlv7YefkUizV6fX\n52XLp9sZ1mhGw+BRAxicza7xzeiolUNyoVimaJmGLNhcKlIb5eFfzuHd594jEUvQd0gfZt10Hgd9\n4cD6cyqG92PT2i00HqsYj8WbXRizNfFYgnVLPyWRUEaMH5q2pzadTl3Yso26yvO3rsSCzaX+fNOD\nrJz/cf2k+c2fbOXOa/7C9X+9mkEjnRkRp196Esvf+ajBFC1/0M9hJ4xtdQpYOpWL1/HHa/9W3zsr\nIlz+319l3JSDM75GodfewAKuGFjngQvt2LSzQajViUVj/CtliaNhY4fwzV9cTJ9BvfD6vPiDPo6e\ncSRf+4+vtPmeNftqueOae9i3q4raqjC1VWFq9tXyp+vvZ/e2vW26VqEODaljHQzul1WNTURuA84G\nIkAlcJmqNr8tuGm3ql3VbP50K30G9Wbbhp3NLnG0ce3mBsfGHXswtz55AzX7ap312jJsOjb2wWvL\nINF0+p1qgoXPL2baRce1+ZqFXHsD62Bws2yboi8BN6pqTER+AdwI3JB9sUydREJ5+JdP8s6TC5ww\ni8YYd8zBRCNpBuv6vIw4tOky4iLSZMnytqraXU08Hm9yPBqOsW9XVbuvW8jP3upYwLlPVk1RVX1R\nVet+wuYB2S8zYRr41/1vMHfOIqLhGDX7aomFYyyf+zH9h/YhENq/tLiIs7/otIumdkg5Dj5qFB5P\n078uwZIAYyZnPyi0UKdlpSr2GQzFJJfP2C4Hnsvh9QzwyoNvNVmtNxqOsnX9dqZ/6xR6VfQkWBpk\n3LFjuP6v36H3gF4dUo4hBw3kyFMmECzZH6bBkgCjJ45s0BObjUJ/9gZdY4pWMWi1KSoiLwPpBjf9\nSFWfTJ7zIyAGPNDCda4ArgAYNHBwuwrbFdXsrUl7PB6Nc+L5Uzj14hM6rSwX/8eXOXTqGObOWUA8\nnuDo6V9g4ukTGkzPyoXUZ2+F2DQF60EtdK0Gm6q2uIaziFwKzACmaQuLu6nqbGA2OOuxta2YXcfK\n+at45k8vs239DoaMGcTAURV8smx9k/P6D+uHP+jvlDJtW7+Df/zySVbOW4XP7+Xo6Udy7r+flXZQ\ncK4U8qyFVI0DzsKtMGTbK3oGcD1wgqpW56ZIXdeilz/kvp/8o37c2Z639+AL+AgE/cRjceLxBB6P\n4Av4ueCGczulTFW7q/nFJb+jem8NmlDisTjvzFnIhtWb+MHd/9bh93dD7Q2sg6HQZNsregcQBF5K\nNkfmqeqVWZfK5eLxBKsWVrJ7616Gjx9KxbDWl3NWVR751ZwGg2lVnV7H4YcOZfCogaxb9hmDRlZw\n2qUnMnjUwI78CvXmPrXQmW+aMtQjFomx/qPP+WT5+ibTsjqCW2pvYFO0CkVWwaaqo3JVkGKxfeNO\nfnPFXVTtqQGFRDzOEdMO42s/OR+Pp/lnUTX7atm3M/2wiY1rt3D9vVd3VJFb9OnKDUTSraIrwqY1\nWzol2Oq4rvZmAZc3NvMgx+6+4T52bd5NuCpMuDpMNBxj8StLefvxd1v8XLAkiK+ZJYLK+/boiKJm\nZOhBgwike5anSsXwzt9YxA3DQupYD2r+WLAl7d62l7cem88bj8xl5+bd7brGzs272VC5mUSjEfqR\n2ghvPDK3xc96fR6O+8oxDcamAQRCAc64/OR2lScXppxzFL6gr0HPpy/gZdDIAQwbl59hi24YFpLK\npmh1Pgs2YN7Ti7hl5i945NdP8dhvnuEn5/0vr/3j7TZfJxKO4pH0f6TR2qZ7CjQ286ozOPaLk/AH\n/QRLAoTKgsy48lSOnn5km8uSK2XlpVz3l6s46KiReLwe/EEfk844gqvv+HrOh3m0lZvCDYp7k5lM\niMgZIvKRiKwWkR+mef9CEVkiIh+KyDsiMiHlvXXJ44tFZGFr9+ryq3vs2rqHv//PY03mXT5++3OM\nnXIw/Yf2zfha/Yb0oaRHSZMBtb6AjyNPm9DMp/bz+jx85dqzOeeq09m3s4qe/bq3e25nLlUM68c1\nd36zfvmjfAdaKjd1LNTpij2oIuIFfg+cCqwHFojIHFVdnnLaWpwRFjtF5Eyc4WFHp7x/kqpuy+R+\nXb7G9sFryxCa/qAm4gnee/nDNl3L4xEuu/UCgiUBvH7neVmwJECfgb049eLMd3MPlgToM6hXQYRa\nKhEpqFBL5bbaG3S5KVqTgNWqukZVI8BDwMzUE1T1HVXdmXyZ1RTNwvrJyYNEPIGSbtUKJZ5mV6jW\nHDTxQH788Pd5+4kFbP98BwdNHMVRp0/otMG0XVmhrxaSTiH3oEYjsbb8Q9G3URNxdnJQfp3BwGcp\nr9fTsDbW2NdpOEVTgZdFJA78sdG1m+jywTb+uEN44ndNd6X3+bwcftKh7bpm7wG9OPvK07ItmmkH\nN6wWkk4hBlzQ56v/88zAtlztBC8iJ+EEW+qKDlNVdYOI9McZN7tSVZvdKbzLN0X7Du7N9G+dhj/o\nx+v14PF6CIT8nDTrWAaP7pxBsCb33Ng0haLuQd0ADE15PSR5rAEROQy4G5ipqtvrjqvqhuR/twCP\n4zRtm9Xla2wAp33tBMZPHcOil5aQiCc4Ytp4hh6c3SbE2UoklBXzV7Fy/iq69ypj0llHUt4vf+PZ\n3MiNTdM6RdjBsAAYLSIjcALtAuCrqSeIyAHAY8DFqvpxyvEywKOqe5O/Pw24taWbWbAlDTywghnf\nOjXfxQCczVbuuOYe1n34GeHqMP6gj2f/9C8uufV8tm3Yycr5q+gzqBcnnD+lQ3aVKiZubZrWKZYp\nWsnFaK8GXgC8wD2qukxErky+fxdwC9AHuDPZSRVLNm8rgMeTx3zAg6r6fJrb1LNgK0Dzn3mPdUs+\nrd+ZvW4oyt03PIA/4HPGy3k9vPvse1z2s68y4YSxgLNz1dqln+IP+Bh+6AEtTuHqatwyHSudQnz+\n1h6q+izwbKNjd6X8/hvAN9J8bg3Q+nipFBZsBWjeM4vqQy2VqtbP20zEE0TiCR746SMcOvVmFr+6\njPtvfRgRp/soVBLg27+9nKFj8tukLiRuHPOWqlgCrjN0+c6DQlQ3Bi4T0XCUlfNXcd9//INwtbN7\nVLgqzO5te7n9qj+l3RshG9FwlEUvfsC/7n+TysXrmuxZ6gZu7VioU8QdDDljNbYCNHXmpAZN0ZbE\n4wmWvrGceJod3+OxOMvnflzfVM3W5nVb+dU3/0A0EiMWieHzeRl+6AFcdftlBTeYuDVubprWKcIO\nhpyxGlsBOuKUwzh82ngCQT++gI9gWZBA0I8/2DA8vF4Pw8YOIRqJEY813UFKFar35G79z7t/eD9V\nu6oIV4WJR+OEayKsXfIJrzz4Vs7u0ZnqJtO7YaWQlnSxGQwZsWArQB6PcMlPzue6v17Fef8+na/e\neC4/f+nHHP+lyfgCPkJlQYIlASpG9OcbP7+IcVMPSbtMdzwe56CJI3NSpl1bdrP5s200bnlGwlHm\nzml1TnJBc3vTFKx52pi72g9dzOBRAxuslPul75/NKV87kU+Wr6e8fw+GHjwIEeGw48cy9OBBfLpi\nQ/0E/GBJgKlfmkyfgbnZtSqRSKSZUZt8L00z2G2KoWkK+2tvXZ0Fm8v07Nudw44/pMExr8/Dv9/5\nDd599n0WvLDYCbVzj2bcsQfn7L69KsrpPaCczZ80XFzBH/Qx6awjcnaffHJ7r6nZz4KtSPj8PqbM\nPIopM4/qkOuLCJf/7EJ+8627nKEmtVGCpUEqDujLKZ24BWBnKJbaW1dmwWYyNnTMIP7r6RtZ8Pxi\ndm3exfDxwzh06hi83uJ7VGvh5m4WbKZNSruXcMJXjsl3MTqFNU3dq/j+qTUmx4qh17SrsWAzJgMW\nbu5iwWZMhty09V9XZ8FmTBu4beu/rsqCzZh2sHArbBZsxrSThVvhsuEexmTBhoQUJquxGZMDVnsr\nLBZsxuSIhVvhsGAzJodsSEhhsGAzJsdsSEj+WbAZ00Es3PLHgs2YDmThlh823MOYDmZDQjpfTmps\nIvIDEVER6ZuL6xlTjKz21nmyDjYRGQqcBnyafXGMKW4Wbp0jFzW23wDX42xAboxphYVbx8vqGZuI\nzAQ2qOoHIs3tYVR/7hXAFQCDBg7O5rbGuJ49d+tYrdbYRORlEVma5tdM4CbglkxupKqzVXWiqk7s\n3btPtuU2pihY7a1jtFpjU9VT0h0XkfHACKCutjYEeE9EJqnqppyW0pgiZhvH5F67n7Gp6oeq2l9V\nh6vqcGA9cKSFmjFtZzW33LJxbMYUCHvuljs5m3mQrLlta/1MY0xLrPaWPZtSZUwBsnDLjgWbMQXK\nwq397BmbMQXMnru1j9XYjHEBq721jQWbMS5h4ZY5CzZjXMTCLTP2jM0Yl7Hnbq2zGpsxLmW1t+ZZ\nsBnjYhZu6VmwGeNyFm5N2TM2Y4pA4+duXZ3V2IwpInUB19VZsBlTZCzcLNiMMZ1ERM4QkY9EZLWI\n/DDN+yIityffXyIiR2b62cYs2IwxHU5EvMDvgTOBscAsERnb6LQzgdHJX1cAf2jDZxuwYDPGdIZJ\nwGpVXaOqEeAhYGajc2YCf1PHPKBcRAZm+NkG8tIrunTZkm0HjRvySQdcui/glsUu3VRWcFd53VRW\n6JjyDsv2AkuXLXnhoHFDMt0EPSQiC1Nez1bV2SmvBwOfpbxeDxzd6Brpzhmc4WcbyEuwqWq/jriu\niCxU1Ykdce1cc1NZwV3ldVNZoXDLq6pn5LsM7WXj2IwxnWEDMDTl9ZDksUzO8Wfw2QbsGZsxpjMs\nAEaLyAgRCQAXAHManTMH+Fqyd3QysFtVN2b42QaKrcY2u/VTCoabygruKq+bygruK2+bqWpMRK4G\nXgC8wD2qukxErky+fxfwLHAWsBqoBi5r6bMt3U9UtcO+jDHG5IM1RY0xRceCzRhTdIoy2ETkByKi\nIpLpGJy8EJHbRGRlcvrI4yJSnu8yNdbWqSz5JCJDReRVEVkuIstE5Jp8l6k1IuIVkfdF5Ol8l6WY\nFF2wichQ4DTg03yXJQMvAYeq6mHAx8CNeS5PA+2ZypJnMeAHqjoWmAxcVeDlBbgGWJHvQhSbogs2\n4DfA9UDB94qo6ouqGku+nIczPqeQtHkqSz6p6kZVfS/5+704gTE4v6VqnogMAaYDd+e7LMWmqIJN\nRGYCG1T1g3yXpR0uB57LdyEaaW6KS8ETkeHAEcD8/JakRb/F+Uc4ke+CFBvXjWMTkZeBAWne+hFw\nE04ztGC0VF5VfTJ5zo9wmlEPdGbZipWIdAMeBb6rqnvyXZ50RGQGsEVVF4nIifkuT7FxXbCp6inp\njovIeGAE8IGIgNOse09EJqnqpk4sYgPNlbeOiFwKzACmaeENKsxkGkxBERE/Tqg9oKqP5bs8LTgW\nOEdEzgJCQA8RuV9VL8pzuYpC0Q7QFZF1wERVLdhVHkTkDODXwAmqujXf5WlMRHw4nRrTcAJtAfDV\n1kZ954s4/6L9Fdihqt/Nd3kylayxXauqM/JdlmJRVM/YXOgOoDvwkogsFpG78l2gVMmOjbqpLCuA\nfxZqqCUdC1wMnJz881ycrBGZLqZoa2zGmK7LamzGmKJjwWaMKToWbMaYomPBZowpOhZsxpiiY8Fm\njCk6FmzGmKLz/wHSTqQT3IY4SAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f67c0be2e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#使用Laplace近似求解bayesian logistic 回归的参数\n",
    "x_train, y_train = create_toy_data()\n",
    "x1_test, x2_test = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))\n",
    "x_test = np.array([x1_test, x2_test]).reshape(2, -1).T\n",
    "\n",
    "feature = PolynomialFeatures(degree=1)\n",
    "X_train = feature.transform(x_train)\n",
    "X_test = feature.transform(x_test)\n",
    "\n",
    "model = BayesianLogisticRegressor(alpha=1.)\n",
    "model.fit(X_train, y_train, max_iter=1000)\n",
    "y = model.proba(X_test)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)\n",
    "plt.contourf(x1_test, x2_test, y.reshape(100, 100), np.linspace(0, 1, 5), alpha=0.2)\n",
    "plt.colorbar()\n",
    "plt.xlim(-5, 5)\n",
    "plt.ylim(-5, 5)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.5 贝叶斯logistic回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算后验概率分布需要对先验概率分布于似然函数的乘积进⾏归⼀化，⽽似然函数本⾝由⼀系列logistic sigmoid函数的乘积组成，每个数据点都有⼀个logistic sigmoid函数。对于预测分布的计算类似地也是⽆法处理的。\n",
    "这⾥我们考虑使⽤拉普拉斯近似来处理贝叶斯logistic回归的问题\n",
    "选择高斯先验：$$p(w)=N(w|m_0,s_0)$$\n",
    "\n",
    "其中m0和S0是固定的超参数。 w的后验概率分布为\n",
    "$$p(w|t)\\propto p(w)p(t|w)$$\n",
    "似然函数为：$$ p(t|w)=\\prod^N_{n=1}y_n^{t_n}{1-y_n}^{1-t_n}$$\n",
    "对数似然：\n",
    "$$lnp(w|t)=\\frac{1}{2}(w-m_0)^T\\mathbf s_0^-1(w-m_0)+\\sum_{n=1}^N{t_ny_n+(1-t_n)ln(1-y_n)}+\\text {常数}$$\n",
    "其中$y_n=\\sigma(w^T\\phi_n)$,最大化后验概率得到$w_{MAP}$，其定义了高斯分布的均值，由laplace近似可知：\n",
    "$$\\mathbf S_N^{-1}=-\\nabla\\nabla lnp(w|t)=S_0^-1+\\sum^N_{n=1}y_n(1-y_n)\\phi_n\\phi_n^T$$\n",
    "于是后验概率分布的高斯近似为:$$q(w)=N(w|w_{MAP},\\mathbf S_N)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.5.1 预测分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对一个新的$\\phi(x)$，类别$C_1$的预测分布可通过对$p(w|t)$积分得到\n",
    "$$P(C_1|\\phi,t)=\\int p(C_1|\\phi,w)p(w|t)dw\\simeq \\int \\sigma(w^T\\phi)q(w)dw$$\n",
    "然后完成marginalization,也就是计算logistic函数与高斯分布的卷积。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "引入狄拉克函数:$$\\delta(x)=\\begin{cases}+\\infty,&x=0\\\\0,&otherwise\\end{cases}$$\n",
    "导数为：\n",
    "$$\\int_{-\\infty}^{\\infty}\\delta(x)dx=1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "于是有$$\\sigma(w^T\\phi)=\\int\\delta(a-w^T\\phi)\\sigma(a)da$$\n",
    "令$$p(a)=\\int \\delta(a-w^T\\phi)q(w)dw$$\n",
    "预测分布就可以写成:\n",
    "$$\\int\\sigma(w^T\\phi)q(w)dw=\\int\\sigma(a)p(a)da$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意到Delta函数给w施加了⼀个线性限制，因此在所有与ϕ正交的⽅向\n",
    "上积分，就得到了联合概率分布q(w)的边缘分布。由于q(w)是⾼斯分布，因此根据2.3.2节，我\n",
    "们知道边缘概率分布也是⾼斯分布。我们可以通过计算各阶矩然后交换a和w的积分顺序的⽅式\n",
    "计算均值和协⽅差，即:\n",
    "\n",
    "$$\\mu_a=\\mathbb E[a]=\\int p(a)ada=\\int q(w)w^T\\phi dw=w^T_{MAP}\\phi$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{align}\\sigma^2_a&=\\int p(a){a^2-\\mathbb E[a]^2}da\\\\&=\\int q(w)\\{(w^T\\phi)^2-(m^T_N-\\phi)^2\\}dw=\\phi^TS_N\\phi\\\\ \\end{align}$$\n",
    "\n",
    "预测分布就变为:\n",
    "$$p(\\mathcal C_1|t)=\\int \\sigma(a)p(a)da=\\int \\sigma(a)\\mathcal N(a|\\mu_a,\\sigma_a^2)da$$\n",
    "关于a的积分表⽰⼀个⾼斯分布和⼀个logistic sigmoid函数的卷积，不能够解析地求值，但是可用逆probit函数Φ(a)的⾼度相似性来获得⼀个较好的近似"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
